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Abstract 

This  report  presents  and  discusses  a  number  of  recent  developments  in  the  steady-state 
thermal  analysis  of  multiple-layer  structures.  These  include:  1)  analytical  evaluation  of 
line  and  area  average  temperatures  and  2)  a  recursion  relation  technique  for  calculating 
the  steady-state  surface  temperature  of  a  multilayer  structure  with  an  arbitrary  number 
of  layers.  The  application  of  the  analytic  averaging  to  the  TXYZ  code  is  incorporated  in 
the  updated  code,  TXYZ30,  while  the  multilayer  recursion  relation  solution  along  with 
analytic  averaging  are  included  in  the  Thermal  MultiLayer  code,  TML.  Both  of  these  are 
contained  in  the  HOTPAC  software  package. 

The  first  part  of  this  report  contains  a  discussion  of  the  general  elements  of  the  multiple- 
layer  thermal  model.  This  is  presented  in  some  detail  for  the  case  of  the  Kokkas  three-layer 
problem  and  the  associated  TXYZ  code.  The  previous  update  of  the  code,  TXYZ20,  is  also 
discussed.  This  incorporates  more  flexible  handling  of  input  data,  assignment  of  positive 
or  negative  noninteger  weights  to  the  various  heat  sources  or  heat  sinks,  and  improved 
evaluation  of  limiting  forms  in  the  code. 

The  second  part  of  the  report  presents  the  analytical  evaluation  of  line  and  area  averages 
which  can  be  calculated  directly.  The  analytical  calculation  of  the  area  average  temper- 
ature should  provide  for  a  more  direct  and  convenient  connection  to  the  experimentally 
measured  temperature  which  tends  to  average  over  the  measurement  area.  The  line  and 
area  averaging  are  incorporated  into  the  TXYZ  code  to  produce  the  TXYZ30  update. 

The  third  portion  contains  a  detailed  discussion  of  the  calculation  of  the  surface  temper- 
ature of  a  multilayer  structure  with  an  arbitrary  number  of  layers.  This  is  based  on  a 
recursion  relation  technique  previously  employed  in  electrical  spreading  resistance  analysis 
to  determine  the  surface  potential  from  the  multilayer  electrical  Laplace  equation.  The 
line  and  area  averaging  techniques  can  be  applied  directly  to  the  recursion  relation  solution 
and  are  included  in  the  Thermal  MultiLayer,  TML,  code. 

The  appendices  contain  the  listing  of  the  annotated,  internally  documented  FORTRAN 
source  codes  for  TXYZ30  and  TML.  These  codes  as  well  as  several  sample  input  and  output 
data  files  are  available  in  ASCII  format  on  DOS-formatted  floppy  disks.  The  sample  input 
and  output  data  files  are  included  so  that  the  user  can  check  the  programs  for  proper 
operation  as  well  as  become  familiar  with  the  setup  and  use  of  the  codes.  Users  of  the 
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previous  versions  of  the  TXYZ  code  should  find  the  TXYZ30  and  TML  codes  easy  to  use 
and  should  benefit  from  the  wider  range  of  problems  which  they  can  be  used  to  address. 

Key  words:  average  temperature;  FORTRAN  programs;  Fourier  analysis;  Laplace  equa- 
tion; multilayer  model;  semiconductor  devices;  semiconductor  materials;  steady-state  heat 
flow;  thermal  conductivity. 
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PART  1.  GENERAL  ELEMENTS  OF  THERMAL  MODEL 

1.1  -  INTRODUCTION 

The  operation  of  a  semiconductor  device  relies  on  the  passage  of  current  through  portions 
of  the  structure.  This  is  accompanied  by  power  dissipation  and  heating  which  gives  rise 
to  a  temperature  distribution.  In  time,  thermal  stresses  may  arise  and  lead  to  device 
degradation  and  possible  failure.  In  addition  to  the  effects  of  thermal  stresses  on  device 
reliability,  the  modeling  of  the  steady-state  thermal  response  of  a  system  is  also  useful 
since  the  temperature  distribution  may  have  an  effect  upon  the  local  mechanical,  thermal, 
electrical,  or  optical  properties  of  the  system.  Consequently,  an  accurate  physical  model 
of  the  temperature  distribution  under  the  power  condition  of  actual  operation  is  of  great 
importance.  It  is  also  possible  to  ascertain  the  relative  effects  of  composition  (thermal 
conductivity)  and  geometry  (layer  thickness).  By  carefully  optimizing  composition  and 
geometry,  it  may  be  possible  to  minimize  the  thermal  stress  and  hence  ensure  optimal 
device  lifetime. 

The  physical  and  mathematical  model  used  here  is  taken  from  the  work  of  Kokkas  which 
has  been  used  previously  to  construct  the  TXYZ  programs.  The  TXYZ  computer  program 
has  been  used  for  a  number  of  years  for  the  thermal  analysis  of  semiconductor  devices  and 
packages.  This  program  makes  use  of  the  closed  form,  Fourier  series  solution  of  the  steady- 
state  heat  flow  equation  for  the  general  case  of  a  rectangular  three-layer  structure  with 
multiple  heat  sources  on  the  top  surface.  TXYZ  provides  for  the  calculation  of  the  tem- 
perature at  any  set  of  points  in  this  structure  and  has  proven  useful  for  the  determination 
of  the  steady-state  temperature  distribution  of  semiconductor  chips  and  packages. 

The  purpose  of  this  work  is  to  report  on  several  recent  advances  which  have  been  made 
in  the  thermal  analysis.  These  include  the  analytical  evaluation  of  line  and  area  averages 
as  well  as  the  recursion  relation  solution  of  the  steady-state  surface  temperature  of  a 
mulitlayer  structure  with  an  arbitrary  number  of  layers. 

1.2  -  STEADY-STATE  HEAT  FLOW:  SINGLE  RECTANGULAR  LAYER 

Consider  a  material  of  uniform  thermal  conductivity,  K\ ,  in  the  form  of  a  rectangular  box  of 
lateral  dimensions  Lx,Ly,  and  thickness  L\ .  The  problem  is  to  determine  the  temperature, 
T(z,  y,  z),  inside  the  material.  The  temperature  is  assumed  to  satisfy  the  steady-state  heat 
flow  equation  [1] 

V2T(z,  *,,.,)  =0.  (1) 

As  this  equation  is  second  order  in  the  three  coordinates,  there  are  six  boundary  conditions. 
Four  of  these  will  be  provided  by  the  lateral  boundary  conditions.  In  the  present  problem, 
all  four  of  the  lateral  boundary  conditions  are  provided  by  the  assumption  that  the  lateral 
surfaces  are  adiabatic;  i.e.,  there  is  no  heat  flow  out  of  the  lateral  boundaries  of  the  material; 
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dT(x,y,z) 
dx 


x=0,L, 


dT(x,y,z) 
dy 


=  0. 


(2) 


y=0,Ll 


The  remaining  two  boundary  conditions  will  be  provided  by  the  vertical  boundary  condi- 
tions (in  z).  These  vertical  boundary  conditions  will  not  be  specified  at  the  present  time 
as  the  intent  of  the  present  section  is  to  obtain  a  general  solution  of  the  one-layer  problem 
where  only  the  lateral  boundary  conditions  are  specified.  As  the  above  equation  is  formu- 
lated in  Cartesian  coordinates,  it  is  convenient  to  use  Fourier  analysis  techniques  to  solve 
the  x  and  y  portion  of  the  equation.  The  Fourier  transform  with  respect  to  the  variables 
x  and  y  is  used,  remembering  that  the  geometry  is  constrained  to  Q,LX,  and  0,Ly.  This 
is  defined  as  [2] 


UxJy,z)=  l       I     T(x,y,z)exp(-2iri(xfx  +  yfy))dxdy, 
Jo  Jo 


(3) 


where  fx,fy  are  the  Fourier  transform  variables  which  are  conjugate  to  the  variables  x,y. 
The  inverse  Fourier  transform  is  defined  as 


/•  +  oo  />+oo 
/       T{fx,fy,z)exp(27ri(xfx  +  yfy))dfxdfy. 
-oo    J  —oo 


(4) 


The  requirement  that  there  is  no  heat  flow  out  of  the  sides  of  the  structure,  i.e.,  dT(x,  y,  z)/dx 
and  dT(x,y,  z)/dy  are  zero  when  x  and  y  are  either  equal  to  zero  or  to  Lx  or  Ly,  respec- 
tively, leads  to  a  consideration  of  the  expression 


dT(x,y,z) 
dx 


/        /       r(fx,fy, z)exp(2Kiyfy)2Trfxl -sin(27rfxx) +  icos(27rfxx)>dfxdfy,  (5) 

J —  oo    J —oo 


and  a  similar  expression  for  dT(x,y,  z)/dy.  If  this  is  to  be  zero  at  the  origin,  this  would 
require  that  the  cosine  term  be  removed.  Next,  consider  the  resulting  expression  at  the 
other  lateral  boundary,  i.e.,  at  x  =  Lx  where  it  is  supposed  to  be  zero.  The  only  way  that 
this  could  be  the  case  is  if  the  argument  of  the  sine  function  is  an  integer  times  7r,  or  that 
the  Fourier  transform  variable  is  of  the  form 
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The  same  argument  applies  to  the  y-dependent  portion.  The  Fourier  representation  of  the 
temperature  with  the  above  lateral  boundary  conditions  may  then  be  written  as 


/+OO  /'+00 
/       T(n,m,z)cos(n7cx/Lx)cos(mTry/Ly)dfxdfy.  (7) 
-oo    J  —  oo 

The  Fourier  transform  equation  is  now  written  as 

r(n,m,2)=  /       /      T(x,  y,  z)  cos(mrx/Lx)  cos(rmry  /  Ly)dxdy.  (8) 
Jo  Jo 

As  the  system  is  of  finite  size,  it  is  convenient  to  write  the  integral  in  eq  (7)  as  a  sum  over 
the  Fourier  cosine  terms  which  fit  into  the  rectangular  geometry.  Further,  as  the  cosine 
function  is  symmetric  around  the  origin,  the  sums  may  be  written  over  only  the  positive 
values  of  the  indices.  It  is  important  to  keep  in  mind  that  the  terms  corresponding  to 
m,n=0  do  not  have  the  factor  of  2  coming  from  the  symmetry  of  the  cos  function.  In 
addition,  the  differentials  may  be  written  as 

»  n       n  + 1       n  1 

dfx  =  A—  =  —  —  =  — .  (9) 

±jx  Lix  JUX  JUX 

The  Fourier  representation  of  the  temperature  may  be  written  as 


7Y         x  _  V"  V"  4T(n,m,z)cos(mrx/Lx)cos(miTy/Ly) 

hho  (*»o  +  l)(«mo  +  l)LxLy  '  (1U) 

where  8nn>  is  the  Kronecker  delta  and  is  equal  to  unity  if  n  —  n'  and  zero  otherwise.  By 
substituting  eq  (10)  into  eq  (1)  and  using 

„2      (  d2       d2       d2  \ 

and 

d2 

— — -  cos(nirx / Lx)  =  —(nw/Lx)2  cos(nirx/ Lx),  (12) 
ox1 

and  the  same  relation  for  the  y-dependence,  it  is  straightforward  to  show  that 
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4  cos(  mr  x  /  Lx)  cosimiry/ Lv)  {    .      ._  ,2     .      ,_  d2  1    .  . 

n=0  m=0         v  /  v  /  »  v 

(13) 


As  the  sum  is  zero  for  arbitrary  values  of  the  variables  x  and  y  (and  the  cosine  terms  are, 
in  general,  nonzero),  then  a  necessary  and  sufficient  condition  that  eq  (13)  is  satisfied  is 
that 


{-(nn/Lx)2  -  (rmr/Ly)2  +  J^}r(n,m,z)  =  0.  (14) 

This  differential  equation  may  be  solved  analytically  using  elementary  methods.  If  the 
variable,  7,  is  defined  as 

H(D2+(^)T-  (i5) 

eq  (14)  may  be  rewritten  as 

d2 

r(n,  m,  z)  —  72t(ti,  m,  z)  =  0.  (16) 


<9z2 

The  solution  of  this  equation  is 

r(n,  m,  z)  =  a  cosh(7^)  +  /?  sinh(72),  (17) 

where  the  coefficients  a  and  /?,  which  may  be  functions  of  7,  are  determined  from  the  two 
2-dependent  boundary  conditions. 

The  above  equation  is  the  general  solution  for  the  z-dependent  Fourier  expansion  coeffi- 
cients for  a  single  rectangular  layer.  It  provides  a  convenient  basis  set  for  the  discussion 
of  the  problem  of  a  rectangular  multilayer  structure  where  all  of  the  layers  have  the  same 
lateral  dimensions.  For  the  multilayer  case,  the  solution  in  each  of  the  layers  can  be  ex- 
pressed in  the  form  of  the  above  equation  where  the  coefficients  are  to  be  determined  from 
the  two  z-dependent  boundary  conditions  appropriate  to  each  of  the  layers.  This  is  used 
in  the  next  section  where  the  multilayer  problem  is  discussed  and  then  specialized  to  the 
three-layer  case. 
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Figure  1.  Geometry  of  the  multilayer  structure  in  which  the  steady-state  temperature  is 
calculated. 
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1.3  -  STEADY-STATE  HEAT  FLOW:  MULTILAYER  RECTANGULAR  STRUCTURE 


The  basic  problem  considered  in  this  section  is  the  general  approach  to  the  calculation  of 
the  temperature  in  a  multilayer  structure  depicted  in  figure  1 .  The  layers  are  characterized 
in  terms  of  the  thermal  conductivities  and  thicknesses  Kl:L{(i  =  1,2,  .  .  .,iV),  where  the 
numbering  begins  at  the  bottom  layer  and  proceeds  up  to  the  top  or  N-th  layer.  Then, 
the  bottom  layer  will  be  layer  number  1  with  temperature  T\(x,y,  z).  The  next  layer  up 
will  be  layer  number  2  with  temperature  T2(x,  y,  z).  This  continues  to  the  top  layer  which 
will  be  layer  N  with  temperature  Tw(x,y,  z). 

The  equations  are  discussed  in  general  for  the  N-layer  problem  and  then  specialized  to  the 
three-layer  problem  as  originally  described  by  Kokkas  [3,4]. 

This  will  allow  for  the  discussion  of  the  Kokkas  three-layer  solution  and  the  TXYZ  code 
while  providing  the  necessary  development  for  the  discussion  of  the  recursion  relation  in 
the  third  part  of  the  report. 

The  thicknesses  and  thermal  conductivities  of  these  layers  are  of  considerable  importance 
in  the  dissipation  of  heat  generated  by  the  power  sources  on  the  surface  of  the  top  layer. 
These  power  sources  are  typically  the  regions  at  or  near  the  surface  of  the  device  where 
currents  are  passed  into  the  device  during  normal  operation.  Consequently,  the  generation 
of  heat  in  the  device  is  one  of  the  unavoidable  side  effects  of  device  operation. 

The  mathematical  formulation  of  this  problem  is  based  upon  the  following  set  of  assump- 
tions: 

1)  the  lateral  dimensions  of  all  the  layers  in  the  structure  are  equal  while  the  thicknesses 
may  be  different; 

2)  each  layer  is  of  uniform,  isotropic,  temperature-independent  thermal  conductivity; 

3)  there  is  no  heat  loss  from  the  lateral  surfaces  due  to  either  radiation  or  convection 
(adiabatic  surfaces)  -  heat  flow  in  the  structure  takes  place  by  conduction; 

4)  there  is  no  input  power  density  inside  the  structure  -  heat  is  generated  only  on  the 
top  surface; 

5)  there  are  no  heat  losses  due  to  interconnections  to  the  top  layer;  and 

6)  the  heat  sink,  which  is  in  contact  with  the  bottom  layer,  is  ideal  and  has  a  temper- 
ature equal  to  ambient. 

The  temperature  in  each  of  the  layers  is  assumed  to  satisfy  the  steady-state  heat  flow 
equation, 

v2r(^,t/,5)  =  o.  (is) 


8 


The  general  one-layer  solution  provides  a  convenient  and  useful  basis  for  the  multilayer 
solution.  Then,  the  temperature  in  each  layer  may  be  written  as 


Ti(x,y,z) 

The  Fourier  coefficients  for  the  solution  of  the  steady-state  heat  flow  equation  in  the  i-th 
layer  of  an  iV-layer  structure  may  be  written  as 

r,(n,  ra,  z)  =  T{(jz)  —  oti  cosh(7z)  +     sinh(72),  (i  =  1, iV),  (20) 

where  the  expansion  coefficients,  a;  and  are  determined  from  the  z-dependent  boundary 
conditions. 

Notice  that  the  assumption  that  there  is  no  heat  flow  out  of  the  lateral  boundaries, 


-0,     (i  =  l,2,  .  .  .,i\0,  (21) 

is  already  contained  in  eq  (19)  as  discussed  in  the  previous  single-layer  analysis. 

The  2N  ^-dependent  boundary  conditions  used  to  solve  the  set  of  equations  in  eqs  (19) 
and  (20)  may  be  expressed  as  follows. 

First,  heat  enters  (or  leaves)  the  structure  through  the  portions  of  the  top  layer  where 
power  is  applied.  This  is  expressed  as 


=  P(*,y),  (22) 

2  =  0 

where  kn  is  the  thermal  conductivity  of  the  top  layer,  and  where  P(x,y)  is  the  power 
function.  This  is  expressed  as  P(x,t/)  =  PoU(x,y),  where  Pq  is  the  power  density  (W/cm2), 
and  U(x,y)  describes  the  surface  geometrical  distribution  of  the  power  sources.  Next,  the 
bottom-layer  boundary  condition  is  provided  by  the  requirement  that  the  temperature  is 
continuous  across  the  interface  between  the  bottom  layer  and  the  heat  sink  and  equal  to 
the  ambient  which  is  taken  as  zero, 


oo 


.  4Tj(n,  m,  z)  cos(mrx / Lx)  cos(rmry / Ly) 

07^0  '  («nO+l)(«mfl+lM  ' 


dTt(x,y,z) 


dx 


dTi{x,y,z) 


x=0.L, 


dy 


dTN(x,y,z) 


Ti(x,y,z)  =  Ta  =  0, 
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(23) 


where  the  function  is  evaluated  at  the  interface.  Finally,  the  remaining  2(7V  —  1)  conditions 
are  provided  by  the  requirement  that  the  temperature  and  the  normal  component  of  the 
heat  flow  are  continuous  across  the  internal  interfaces.  These  are  expressed  as 

Ti(x,  y,  z)  =  T?_j  (x,  y,  z),  (24) 


oz  oz 

where  the  functions  and  their  derivatives  are  to  be  evaluated  at  the  interfacial  boundaries. 

For  the  case  of  an  TV-layer  structure,  the  substitution  of  eqs  (19)  and  (20)  into  the  boundary 
conditions  given  by  eqs  (22  -  25)  gives  rise  to  a  set  of  2N  equations  in  2iV  unknowns 
({a;(7)},  {fliij)},  i  —  l,...,iV).  The  analytic  solution  of  this  system  of  equations  requires 
the  use  of  matrix  algebra  involving  Cramer's  rule  [5]  and  the  Laplace  method  [5]  for  the 
evaluation  of  the  resulting  determinants.  Clearly,  this  also  can  become  rather  tedious, 
especially  since  the  expansion  coefficients  are  functions  of  the  continuous  variable,  7. 

For  cases  up  to  i  =  3,  Kokkas  [3,4]  was  able  to  work  out  the  system  of  equations  for  not 
only  the  surface  temperature  but  also  the  temperature  inside  the  three-layer  structure. 
The  Kokkas  three-layer  solution  is  reviewed  in  the  next  section  and  the  N-layer  problem 
is  returned  to  in  Part  3. 

1.4  -  KOKKAS  THREE-LAYER  MODEL 


(25) 


The  Kokkas  three- layer  model  [3,4]  involves  the  solution  of  the  multilayer  equations,  eqs 
(19,  20,  22  -  25),  for  the  case  where  N=3.  This  can  be  carried  out  analytically  as  the 
systems  of  equations  are  still  tractable.  Specializing  the  multilayer  equations  to  the  case 
of  three  layers  yields  the  Fourier  coefficients 


Tj(n,  m,  z)  =  Ti(jz)  =  at  cosh(7.?)  +  fa  sinh(72),  (i  =  1,2, 3). 


(26) 


The  corresponding  boundary  conditions  are 


*3 


dTN(x,y,z) 


z=0 


T3(x,y,z)  =T2(x,y,z)\ 


T2(x,y,z) 


z=-L 


z=-(L3+L2) 


=  Ti(x,y,z) 


z=  —  L3 


z=-(L3  +  L2) 


(27) 
(28) 
(29) 
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«3 


dT3(x,y,z) 
dz 


z=  —  Lz 


«2 


dT2(x,y,z) 
8z 


(30) 


Z=  —  L; 


*2 


dT2(x,y,z) 
dz 


=  «1 


z=-(L3  +  L2) 


dTx{x,y,z) 
dz 


(31) 


z=-(L3  +  L2) 


Ti(x,y,z) 


z—  —  L, 


=  Ta  =  0, 


(32) 


where  all  temperatures  are  measured  relative  to  the  ambient  heat  sink  temperature  and 
Lz  —  L$  +  L2  +  L\. 

It  is  important  to  note  that  the  origin  of  the  depth  scale  is  at  the  surface  of  the  top  layer 
and  that  all  vertical  distances  are  negative. 

The  solution  of  the  system  of  equations  is  by  standard  matrix  algebra  techniques.  However, 
instead  of  using  this  method  explicitly  here,  the  Fourier  coefficients  for  each  of  the  layers  are 
presented  and  are  shown  to  satisfy  the  heat  flow  equation  and  the  appropriate  boundary 
conditions.  In  particular,  the  Fourier  coefficients  in  the  three  layers  are  presented  in 
reference  [4].  Specializing  these  to  the  steady-state  situation,  they  are 


r3(n,m,2)  =  Aj£cosh(7(L3  +  z))  +  Csinh(7(L3  +  z))j, 


(33) 


where 


r2(n,m,z)  =  A^D  cosh(7(L3  +  L2  +  z))  +  E  sinh(7(L3  +  L2  +  z)) },  (34) 


Ti(n,m,  z)  =  Asinh(7(L3  +  L2  +  L\  +  z)), 


(35) 


A  = 


U(n,  m)Po 


k37       (  B  sinh(7L3)  +  C  cosh(7L3) 


}■ 


(36) 


B  =  Dcosh(7L2)  +  £sinh(7L2), 


(37) 


(38) 
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D  =  sinh(7l/i), 


(39) 


E  =  —  cosh(7jLi),  (40) 

«2 


and 

'L>x  r^y 


U(n,m)  =  /       /      U(x,y)cos(nirx/Lx)cos(miry/ Ly)dxdy,  (41) 
Jo  Jo 

is  the  double  Fourier  cosine  transform  of  the  power  density  uniformity  function. 

Now,  it  is  shown  that  the  above  are  the  solutions  of  the  z-dependent  part  of  the  steady-state 
heat  flow  equation  (see  eq  (16)) 

d2 

Ti(n,m,z)  -  j2Ti(n,m,z)  =  0,  (42) 


d 

where  the  subscript  i  takes  on  the  values  of  3,2,1.  This  may  be  easily  shown  to  be  the  case 
as 


d2 

cosh(7(L  +  z))  =  72  cosh(7(£  +  z)) ,  (43) 


2 


dz 


and 

d2 


dz< 


sinh(7(i:  +  z))  =  72  sinh(7(L  +  *)),  (44) 


where  L  is  a  constant  and  is  equal  to  L3,  L3  +  X2,  or  £3  +  £2  +  L\  in  eqs  (33  -  35).  Hence, 
eqs  (33  -  35)  satisfy  the  z-dependent  differential  equation.  Next,  it  is  shown  that  these 
Fourier  coefficients  satisfy  the  appropriate  boundary  conditions.  The  first  of  these  is  that 


drgfu,  m,  z) 
«3- 


dz 


=  P(n,  m)  =  U(n,  m)P0 .  (45) 

2  =  0 


Using  the  Fourier  coefficient  given  by  eq  (33)  for  the  top  layer,  this  may  be  evaluated  as 
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«3 


dr3(n,m,z) 


2  =  0 


ac37^4<  Z?sinh(7.L3)  +  C  cosh(7.L3) 


«37 


U{n,  m)P0 


K37       [  5sinh(7L3)  +  Ccosh(7L3) 


#sinh(7L3)  +  C  cosh(7i:3) 


) 


f7(n,  m)Po 


(46) 


Hence,  the  top-layer  boundary  condition  is  satisfied  by  the  T3(n,m,z).  Next,  consider  the 
bottom-layer  boundary  condition,  i.e., 


Ti(n,m,z) 


z=-(L3  +  L2+Ll) 


(47) 


Making  use  of  eq  (35),  this  may  be  readily  evaluated  as 


ri(n,m,z) 


z=-(L3  +  L2+L1) 


=  A  sinh(7(L3  +  L2  +  £1  -  L3  -  L2  -  £1))  =  0.  (48) 


Hence,  the  last  boundary  condition  is  satisfied.  The  final  set  of  boundary  conditions  to  be 
verified  are  the  ones  which  pertain  to  the  interface  boundary  conditions.  These  are 


T3(n,ra,2)  =r2(n,m,2) 


z=-L 


z=  —  L3 


(49) 


T2(n,m,z) 


z=-(L3+L2) 


Ti(n,m,z) 


z=-(L3  +  L7) 


(50) 


«3 


dT3(n,m,z) 


dz 


z=  —  L3 


K2 


dT2(n,m,  z) 


z—  —  L3 


(51) 


«2 


dT2(n,m,z) 


=  «i 


z=-(L3+L2) 


dr\{n,  m,  z) 


z=-(L3  +  L2) 


(52) 


In  order  to  verify  these  equations,  it  is  simplest  to  calculate  the  right-  and  left-hand  sides  of 
the  equations  and  then  compare  them  directly.  The  left-hand  side  of  the  first  temperature 
continuity  equation  may  be  evaluated  as 
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T3(n,m,z) 


z=  —  L3 


A^Bcosh(j(L3  -  L3))  +  Csinh(7(L3  -  Ls))} 
=  AB  —  a|d  cosh(7L2 )  +  E  sinh(7i:2 ) } • 


(53) 


The  right-hand  side  of  the  equation  may  be  evaluated  as 


r2(n,m,z) 


2=  — Ls 


=  A|JDcosh(7(L3  +  L2  -L3))  +  Esinh(7(Z-3  +  L2  -£3))} 

=  A{D  cosh(7£2 )  +  E sinh(7i:2 ) } •  (54) 


Hence,  the  first  of  the  temperature  continuity  equations  satisfies  the  boundary  condition. 
Next,  consider  the  second  temperature  continuity  equation, 


T2(n,m,z) 


z=-(L3+L2) 


=  TX(n,m,z) 


z=-(L3  +  L2) 


(55) 


The  left-hand  side  of  the  equation  may  be  evaluated  as 


r2(n,m,z)  =  A\  D  cosh(j{L3  + L2- L3- L2))  +  E  sinh(y(L3  + L2- L3- L2))\ 

z—-(L3  +  L2)  {  > 


=  AD  =  Asinh(7Li). 


(56) 


The  right-hand  side  is 


ri(n,  m,  z) 


2=-(L3+L2) 


Asmh(j(L3  +  L2  +      -  L3  -  L2))  =  Asinh^Li).  (57) 


Hence,  the  second  temperature  continuity  equation  is  satisifed.  The  next  equation  to  be 
verified  is  the  first  heat-flow  continuity  equation;  i.e., 


«3 


<9r3(n,m,z) 


3; 


«2 


dT2(n,m,z) 


2=  —  L? 


(58) 


z=  —  L3 


Evaluating  the  left-hand  side  leads  to 
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«3 


5r3(n,m,  z) 


dz 


z=-L: 


=  «37^{^smh(7(JL3  -L3))  +Ccosh(7(L3  -L3))} 
=  KzjAC  =  K2jA^Dsmh(-)L2)  +  £cosh(7£2)}.  (59) 


The  right-hand  side  may  be  evaluated  as 


«2 


dr2(n,m,  2:) 


5; 


=  K27^{^sinh(7L2)  +  £cosh(7L2)} 


(60) 


The  final  equation  to  be  evaluated  is  that  for  the  heat  flow  continuity  between  the  second 
and  third  layers;  i.e., 


«2 


dT2(n,m,z) 


dz 


«1 


dri(n,  m,  2) 


2=-(L3+L2) 


z=-(L3  +  L2) 


The  left-hand  side  of  this  equation  may  be  evaluated  as 


«2- 


<9r2(rc,m,  z) 


dz 


2=_(L3+L2) 


K27A|Dsinh(7(i:3  +L2  -  L3  -  L2))  +  E  cosh(7(£3  +  L2  -  L3  -£2))  j 
=  ac27A£i  =  K17A  cosh(7Z/i ). 


The  right-hand  side  is 


(61) 


(62) 


dri(n,  m,  z) 


dz 


=  /c17Acosh(7(L3+L2+Li-L3-JL2)  =  KijAcosh(jLi).  (63) 


2=_(L3+L2) 


Now  that  it  has  been  shown  that  the  Fourier  coefficients  satisfy  the  steady-state  heat  flow 
problem  and  the  appropriate  boundary  conditions,  there  are  several  points  to  be  considered 
before  getting  into  the  body  of  the  program.  These  include:  (1)  evaluation  of  the  function 
J7(n,m)  for  a  uniform  power  source  of  given  size  and  (2)  simplification  of  the  Fourier 
coefficients  for  .subsequent  numerical  analysis.  The  latter  point  is  necessary  as  the  limits 
of  small  7  and  large  7  may  give  rise  to  overflow  or  underflow  problems  when  the  program 
is  constructed. 


15 


1.5  -  KOKKAS  THREE-LAYER  MODEL:  FORM  OF  THE  FUNCTION  U(n,m) 


The  first  thing  to  be  considered  is  the  form  of  the  function  U(n,m)  for  an  arbitrary  number 
of  heat  sources.  This  function  is  defined  as 


l-LX  pLy 

U(n,m)  =  I       I      U(x,  y)  cos(mtx/Lx)  cos(mny  /  Ly)dxdy.  (64) 
Jo  Jo 

The  analysis  can  most  easily  be  accomplished  in  terms  of  a  single  uniform  heat  source.  The 
case  of  an  arbitrary  number  of  heat  sources  follows  by  summing  the  results  of  each  heat 
source.  Further,  if  any  of  the  heat  sources  are  nonuniform,  their  effects  can  be  constructed 
by  suitably  overlapping  a  number  of  uniform  heat  sources.  In  the  coordinate  system  being 
used,  consider  a  single  heat  source  denoted  by  the  index  i  with  a  corner  at  the  location 
(xi,yi)  and  lengths  along  the  x-  and  y-directions  given  by  (lxlJyi).  Over  the  area  of  the 
heat  source,  U(x,y)  is  assumed  to  be  uniform  and  equal  to  unity.  Away  from  the  area 
of  the  heat  source,  U(x,y)  is  assumed  to  be  equal  to  zero.  Then,  U(x,y)  may  be  viewed 
as  being  a  unit  step  function  over  the  surface  of  the  power  source.  Consequently,  the 
contribution  from  this  single  heat  source  may  be  written  as 


I       I      U(x,  y)  cos(nirx/ Lx)  cos(miry/Ly)dxdy  = 
Jo  Jo 

rxi  +  lxi     ryi  +  lyi 

I  I  cos(nirx / Lx)  cos(mny / Ly)dxdy.  (65) 

The  integrals  can  be  simply  evaluated  to  give  the  result  that 


TT  LxLy     f  .   /mr(xi  +  lx{)\       .  /mrxi\} 

Ui(n,  m)  =  -  -  y  -  <  sin  — ^—   -  sin  — —  > 

v  '    ;      (n7r)(m7r)  \     V        Lx        )         V  Lx  )  J 

J.   (mir(yi  +  lyi)\       .  /rmryl\\ 

XlSml  Ly  )-Sm{-Ly-)j-  ^ 

Similar  expressions  may  be  written  for  each  of  the  heat  sources  and  then  summed  to  give 
the  cumulative  heat  source  effect. 

Before  turning  to  the  small  7  and  large  7  behavior  of  the  Fourier  coefficients,  it  is  important 
to  consider  the  behavior  of  the  function  U(n,m)  for  either  n  =  0,  m  =  0,  or  both.  This 
is  important  in  the  numerical  implementation  of  the  solutions  as  the  program  will  have  to 
calculate  the  double  Fourier  cosine  transform  over  the  range  of  n,m  required  by  the  sum 
in  eq  (19).  Once  the  value  of  n  or  m  is  zero,  there  will  be  problems  with  most  machines  as 
far  as  evaluating  the  seeming  divergence.  This  can  be  circumvented  by  investigating  the 
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behavior  of  the  function  for  n  =  0.  This  may  be  readily  carried  out  by  first  noting  that  the 
function  ?7(n,  m)  (for  the  particular  form  of  U(x,  y))  is  a  product  of  two  terms.  This  may 
be  simply  written  as  U(n,  m)  =  U(n)U(m).  Then,  consider  the  n  (from  the  x  integration) 
contribution  to  the  function  which  is  given  by 


U[n)  =  _|sm^  _  }_sm(__jj.  (67) 


There  is  an  apparent  divergence  or  infinity  if  n  is  simply  set  equal  to  zero.  This  is  the  way 
in  which  a  computer  would  look  at  the  expression.  However,  this  infinity  is  not  real  as  can 
be  seen  by  using  the  expansion  of  the  sine  function  for  small  values  of  the  argument.  In 
particular, 

x3 

sm(x)  =  x  -  —  -\  .  (68) 

Making  use  of  this  expression  for  the  sine  function,  it  is  straightforward  to  show  that 


v     tt<  \      v     L*  f  •   fnir(xi +lxi)\  fnnxt\\ 
hm  U  ( n )  —  urn  —  <  sin    —  sin    >  = 

n_0  n— 0  717T  [        V  Lx  /  V    Lx    J  J 

hstl — lx — )  ~  Kir))  =  lXi-  (69) 

The  same  conclusion  holds  for  the  y-dependent  portion,  i.e.,  U(m).  This  must  be  specially 
coded  to  bypass  any  overflow  problem.  The  specific  coding  of  the  heat  source  may  be 
found  in  the  program  listing  in  the  function  UZERO. 

In  general,  the  function  U(n,m)  is  oscillatory  and  does  not  approach  zero  sufficiently  fast 
for  large  values  of  n  or  m.  In  particular,  eq  (66)  shows  that  Ul{n,m)  — >  0  like  1/nm  as 
m,n  — >  oo.  In  addition,  the  cosine  terms  in  eq  (19),  i.e.,  cos(nnx/ Lx)  cos(m7ry/L2/),  do 
not  approach  a  definite  limit  for  large  values  of  the  argument.  Consequently,  the  product 
U(n,  m)  cos(mvx / Lx)  cos(rmry / Ly)  tends  to  zero  slowly. 

1.6  -  KOKKAS  THREE-LAYER  MODEL:  BEHAVIOR  OF  FOURIER  COEFFICIENTS 
FOR  SMALL  VALUES  OF  THE  ARGUMENT 

As  has  been  seen  in  the  treatment  of  the  function  U (n,  m),  care  must  be  taken  for  the  case 
where  both  n  and  m  are  equal  to  zero  (or  7  =  0).  The  same  considerations  must  be  carried 
out  for  the  Fourier  coefficients  in  the  three  layers.  In  the  solutions  in  the  three  layers,  the 
summation  indices  n,m  appear.  The  summation  over  these  variables  is  of  influence  in 
the  variable  7  according  to  eq  (15).  Also,  the  Fourier  coefficients  contain  the  hyperbolic 
functions  which  depend  upon  7.  For  large  values  of  7,  the  sinh  and  cosh  functions  grow 
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exponentially.  This  can  present  special  numerical  problems  when  the  summation  variables 
approach  the  upper  limits  which  may  be  required  for  the  case  of  very  small  heat  sources. 
Hence,  special  care  must  be  taken  to  study  the  behavior  of  the  Fourier  coefficients  for  small 
7  and  large  7  to  remove  any  potential  numerical  overflow  problems.  Once  this  is  properly 
taken  care  of,  the  Fourier  coefficients  and  the  solutions  will  be  numerically  well  behaved. 

First,  consider  the  small  7  behavior  of  the  Fourier  coefficients.  This  is  done  by  considering 
the  small  7  behavior  of  the  eqs  (33  -  40)  and  the  small  argument  behavior  of  the  hyperbolic 
functions.  In  the  following  discussion  as  well  as  the  discussion  of  the  large  7  behavior,  the 
term  U(n,  ra)Po / k3  will  be  removed  for  convenience.  This  term  will  henceforth  be  included 
explicitly  in  the  sum  in  eq  (10)  as  the  Fourier  coefficients  for  all  three  layers  contain  this 
as  a  common  factor  through  A  (see  eqs  (33  -  36)).  Then,  for  small  7, 


E 


11 

«2 


(70) 


D  «  jLi ,  (71) 


C^M  +  -L  (72) 

K3  I  K2  I 


B^1L1  +  — 7L2,  (73) 

K2 


and,  remembering  that  the  factor  U(n,m)Po  / k$  has  been  included  explicitely  in  eq  (19), 


A«i-3-.  (74) 

7  «1 

Making  use  of  these  expressions  and  the  small  argument  behavior  of  the  hyperbolic  func- 
tions, it  is  straightforward  to  investigate  the  small  7  behavior  of  the  solutions.  In  particular, 

.  .         1  K$  [     _  K\  K2  f    9  ,    _  K\  1     ,  _  ,1 

T3(n,m,z)  w  <^  7L1  +  — 7Z2  +  —WL2LX  +  —  WL3  +  z)  } 

j  K-i  {  k2  k3  I  k2)  J 

w  —  i  Li  +  —L2  +  —  (£3  +  z) 

«1  [  «2  «3 

Then, 
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^3  ^3 

lim  r3(n, m, z)  =  (L3  +  z)  H  1-2  H  £i-  (75) 

7^0  «2  /tx 


Next, 


r2(n,m,z)  w  -  —  I  7^  +  —7(^3  +  L2  +  2)}, 
7Ki  I  «2  J 


Then, 


]im  T2(n,m,z)  =  —L1  + —(Ls  +  L2  +  z).  (76) 
7^0  Kl  K2 


And  finally, 


1  ac3 

Ti(rc,m,z)  «  <[  7(^3  +  L2  +  L2  +  z) 

7  AC! 


Then, 


^3 

limri(n,m,2)  =  — (L3  +  L2  +  Li  +  2).  (77) 

7^0  K\ 

These  special  forms  of  the  Fourier  coefficients  (in  the  limit  as  7  — >  0)  are  necessary  in  the 
code  to  bypass  overflow  problems  for  small  values  of  the  argument. 

1.7  -  KOKKAS  THREE-LAYER  MODEL:  BEHAVIOR  OF  FOURIER  COEFFICIENTS 
FOR  LARGE  VALUES  OF  THE  ARGUMENT 

As  the  Fourier  coefficients  have  been  investigated  for  small  7  and  have  been  shown  to  be 
well-behaved  when  properly  written,  what  remains  is  to  write  these  coefficients  in  a  form 
which  is  amenable  for  investigating  their  large  7  behavior.  As  noted  before,  the  hyperbolic 
functions,  sinh  and  cosh,  grow  exponentially  for  large  values  of  the  argument.  On  the 
other  hand,  the  hyperbolic  tanh  approaches  unity  for  large  values  of  the  argument.  With 
this  in  mind,  let  us  investigate  the  form  of  the  Fourier  coefficients,  written  as  much  as 
possible  in  terms  of  the  tanh,  which  takes  care  of  this  potential  numerical  difficulty.  To 
this  end,  it  is  convenient  to  introduce  the  shorthand  notation  for  the  hyperbolic  functions, 
c(x)  =  cosh(x),  s(x)  =  sinh(x),  and  t(x)  —  tanh(z).  Making  use  of  this  shorthand 
notation,  the  Fourier  coefficients  may  be  written  as 

r3(n,m,z)  =  a[Bc(7(L,  +  z))  +  Cs(7(Ls  +  2))},  (78) 
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r2(n,  m,  *)  =  a{Dc(j(L3  +  L2  +  «))  +  £*(7(L3  +  ^2  +  *))  },  (79) 


where, 


T1(n,m,z)  =  As(j{Li  +  L2  +  Li  +  z),  (80) 


A     7t^(7L3)  +  CC(TL3)i'  (81) 
5  =  Dc(iL2)  +  Es{jL2),  (82) 

C=  —  {^(7X2)4-^(7X2)},  (83) 


D  =  j(7Ii),  (84) 


E='^c{1L1).  (85) 

«2 


As  in  the  investigation  of  the  small  7  behavior  of  the  Fourier  coefficients,  the  factor 
U(n,m)Po I K3  has  been  deleted  from  eq  (81)  for  convenience.  This  factor  may  simply 
be  included  in  the  Fourier  representation  of  the  temperature,  eq  (19),  as  it  is  common  to 
all  three  layers. 

Now,  the  above  equations  (eqs  (78  -  85))  will  be  rewritten  by  making  use  of  the  definition 
of  the  hyperbolic  tanh;  i.e.,  t(x)  —  s(x)/c(x).  First,  consider  the  coefficient  C. 


C  =  —{Ds(7L2)  +  Ec(7L2) 

«3 


}■ 


C  =  ^  s(7Li)s(7L2)  +  ~c(7£i)c(7L2)  I, 

«3  I  «2  J 


c  =  ^c^l^lA  ti-yLiX-YU)  +  - 

«3  I  K2 


(86) 
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Next,  the  coefficient  B  may  be  written  as 

B  =  Dc(jL2)  +  Es(jL2), 


B  =  s{1L1)c{lL2)  +  —0(7^)5(7X2), 

«2 


B  =  c(TL1)c(7L2)|<(7L1)  +  ^t{jL2) 


Then, 


Bs(yL3)  +  Cc(1L3)  = 

c{1L1)c{lL2)s{iLAi{1L1)  +  — *(7X2)}  +  0(7X00(7X2)0(7X3) — (tfrl^frLa)  +  —  j 
-  c(7X1)c(7X2)C(7L3)(<(7iiW7^3)  +  —  t(7L3)t(7L2)  +  ^frXj^X,)  +  —  j.  (88) 

[  «2  ^3  «3  J 

Then,  the  coefficient  A  may  be  written  as 


A  = 

1 


(89) 


70(7X^(7X3)0(7X3)  {t{lLl)t{lLz)  +  S,(7X3)<(7X2)  +  S^WTX,)  +  g 
It  is  convenient  to  define  the  function  ^(7)  as 


0(7)  -  7  T~s  (90) 

|<(7X1)*(7X3)  +  %t(<yL3)t(7L2)  +  g*(7X1)*(7X2)  +  ^} 

which  is  well  behaved  for  all  values  of  7.  Then,  the  coefficient  A  may  be  written  as 

A=     .       "'7>       ...  (91) 

70(7X1)0(7X2)0(7X3) 
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By  making  use  of  the  above  procedure,  it  is  relatively  straightforward  to  show  that  the 
Fourier  coefficients  as  given  by  eqs  (78  -  80)  may  be  written  as 


r3(n,m,z)  = 

n(7)c(7(£.  +  »))  f  +  «,  +  <(        +  z))  M  /  }  +  «.)  1  (92) 

70(7X3)  [  K2  AC3  V  «2/  J 

0(7)^(7(^3+^2+^))  f  v    ,  ,r  ,      ,      vx\  ,Q,v 

r2(n,m,z)  =    <  H7^i)  +  — %(£s  +^2  +  2  )  f,  93 

70(7X3)0(7^2)        t  ^2  J 


and 


0(7)^(7(^  +  ^2+^+2)) 

n(n,m,2)  =     -.  (94) 

7c(7L3)c(7L2)c(7Li) 

In  eqs  (92  -  94),  the  function  0(7)  and  the  terms  inside  the  curly  brackets  are  well  behaved 
for  all  values  of  the  variable  7.  The  sinh  and  cosh  terms  which  remain  may  still  give  rise  to 
numerical  overflow  problems  for  large  values  of  the  argument.  However,  as  both  of  these 
functions  grow  exponentially  for  large  values  of  the  argument  and  they  appear  in  both  the 
numerator  and  the  denominator  of  the  Fourier  coefficients,  there  will  be  cancellation.  This 
cancellation  for  large  values  of  the  argument  will  not  be  worked  out  in  detail  here  but  is 
contained  in  the  FORTRAN  listing  of  the  program  in  the  function  FUNZ. 

These  Fourier  coefficients  are  used  in  the  equation 


„              4U  (n,  m^n,  m,z)cos(mrxLx)cos(m'Ky/Ly) 
Ti(x,y,z)  =  P0  >    >   77 — — T77-  ,1U  T   95) 


for  i  =  1,2,3  to  obtain  the  solutions  in  each  of  the  three  layers.  In  the  above,  the  term 
PoU(n,m)/ K3  has  been  written  out  explicitly  and  is  no  longer  contained  in  the  Fourier 
coefficients.  It  is  important  in  obtaining  the  solution  in  x,y,  z  to  use  the  appropriate  layer 
equation.  This  is  automatically  taken  care  of  in  the  program  as  the  depth  z  is  compared 
with  the  various  thicknesses  and  the  corresponding  layer  Fourier  coefficient  is  used.  The 
user  does  not  have  to  specify  which  layer  is  to  be  used. 

There  are  several  parts  of  the  above  equation  which  can  be  identified.  First,  there  is 
the  z-dependent  part  of  the  solution  in  each  layer,  Tj(n,m,z).  The  two  cosine  functions, 
cos(mrx / Lx)  cos(rmry / Ly),  arise  from  the  x-  and  y-dependent  parts  of  the  steady-state 
equation.  Their  specific  form  arises  from  the  requirement  of  the  lateral  boundary  condi- 
tions. Finally,  the  requirement  that  heat  enters  or  leaves  the  top  surface  only  through  the 
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heating  elements  introduces  the  function  PoU(n,m).  It  should  be  noted  that  Pq  appears 
only  in  front  of  the  sum.  It  is  a  scale  or  multiplicative  factor  which  for  convenience  may 
be  set  equal  to  unity.  Its  effect  is  to  scale  the  temperature  in  a  direct  linear  fashion. 
Performing  the  calculation  for  Pq  =  1  for  a  given  structure  then  provides  the  temperature 
distribution  for  arbitrary  Pq  which  can  then  be  obtained  by  multiplying  by  the  Pq  which 
is  used. 

An  interesting  exercise  left  to  the  reader  is  to  show  when  the  three  thermal  conductivities 
are  equal  that  the  one-layer  solution  is  obtained.  Also,  another  exercise  is  to  show  that 
the  one-layer  solution  is  obtained  when  the  thicknesses  of  the  second  and  third  layers  are 
set  equal  to  zero. 

1.8  -  SPECIAL  CASE  OF  POWER  SOURCE  COVERING  TOP  SURFACE 

As  a  special  case  of  eq  (95),  consider  the  situation  of  a  single  power  source  completely 
covering  the  top  surface;  i.e.,  there  is  a  single  heat  source  with  lateral  dimensions  equal 
to  that  of  the  three-layer  structure.  In  this  particular  case,  it  is  shown  that  the  above 
equation  reduces  to  the  familiar  thermal  resistance  equation.  The  easiest  way  to  proceed 
with  the  analysis  is  to  consider  the  specific  form  of  the  function  U(n,m).  From  eq  (66), 


Ui(n,  m)  = 

(96) 

For  the  particular  case  of  uniform  surface  coverage,  X\  =  y\  =  0,  Ixi  —  Lx,  and  ly\  —  Ly. 
Upon  substituting  these  values  into  the  equation,  the  function  reduces  to 


U\(n,m)  =  *   y     |sin(n7r) sin(m7r)  1.  (97) 

(n7r)(m7r)  ^  J 

This  is  zero  when  the  indices  are  nonzero.  For  the  case  where  both  of  the  indices  are  zero, 
the  use  of  the  expansion  of  the  sine  function  gives  rise  to  the  result  that 


Ui(n,m)  =  LxLySn0Smo.  (98) 

If  this  form  of  the  U(n,  m)  function  is  substituted  into  eq  (95)  for  the  temperature  in  each 
of  the  three  layers  and  the  form  of  the  Fourier  expansion  coefficients  as  7  — >  0  (eqs  (75  - 
77))  is  used,  it  is  readily  shown  that  the  temperatures  in  the  three  layers  may  be  written 
as 

I*3(W)  =  Po{^  +  ^  +  H  (99) 
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L3  +  L2  +  z  Li 
T2{x,y,z)  =  P0\  +  — 

K2  Ki 


(100) 


rp  (  v  p  f  U+L2+LX+z\ 
2l(x,y,z)  =  P0|   —   J 


(101) 


As  Pq  is  the  power  density  per  unit  area,  these  equations  give  rise  to  the  usual  results  of 
the  one- dimensional  calculations  of  the  thermal  resistance. 

1.9  -  EFFECT  OF  UPPER  SUMMATION  LIMITS  ON  TEMPERATURE 

It  is  important  to  keep  in  mind  that  the  Fourier  coefficients  are  functions  of  the  variables  n 
and  m.  As  discussed  in  the  previous  sections,  the  power  density  function,  U(n,m),  tends 
to  zero  very  slowly  for  large  values  of  the  argument.  Also,  the  cosine  terms  in  eq  (95) 
do  not  tend  to  any  limit  as  the  arguments  approach  infinity.  It  is  the  Fourier  coefficients 
which  are  responsible  for  the  convergence  of  the  sum.  This  is  especially  the  case  for  the 
surface  (z  =  0)  value  of  the  temperature.  This  consideration  comes  into  play  when  setting 
the  upper  summation  limit. 

Consider  the  case  of  a  structure  with  lateral  dimension  of  L  and  a  heat  element  of  lateral 
dimension  of  A.  For  the  Fourier  series  to  begin  "seeing"  this  element,  the  cosine  function 
in  the  basis  set  must  have  at  least  one  complete  cycle  in  A.  As  the  element  is  rectangular, 
one  cycle  is  certainly  not  enough,  as  the  cosine  is  a  poor  representation  of  the  rectangle. 
Consequently,  higher  "frequencies"  are  required  to  assure  the  adequacy  of  the  representa- 
tion. The  rule  of  thumb  is  that  the  number  of  terms  should  be  at  least  on  the  order  of 
L/A.  A  stronger  rule  of  thumb  would  require  several  times  L/A.  In  general,  the  stronger 
rule  of  thumb  should  be  applied  to  the  smallest  heat  element  to  achieve  better  accuracy. 

1.10  -  THE  TXYZ  AND  TXYZ20  CODES 

The  above  equations  have  been  coded  in  the  TXYZ  program  [6].  Over  time,  various 
updates  were  included  in  the  TXYZ20  program  [7].  In  the  original  version  of  TXYZ,  all 
heat  sources  had  equal  unit  weights.  The  updated  TXYZ20  program  allows  for  noninteger 
weights  to  be  assigned  to  each  heating  element.  This  weight  may  be  positive  or  negative 
and  hence  represent  the  effects  of  either  a  heat  source  or  a  heat  sink. 

The  reason  for  the  relatively  compact  nature  of  the  original  TXYZ  code  and  its  numerical 
efficiency  is  that  great  care  was  exercised  in  the  investigation  of  the  functions  for  the  small 
and  large  argument  regimes.  For  small  values  of  the  argument,  the  evaluation  of  limiting 
forms  led  to  expressions  which  were  correct  and  free  of  artificial  numerical  singularities. 
For  large  values  of  the  argument,  special  care  was  required.  Many  of  the  functions  involved 
contained  the  hyperbolic  sinh  and  cosh  functions.  These  approach  the  exponential  function 
for  large  values  of  the  argument.  These  provide  the  possibility  for  numerical  infinities  or 
overflows.  Careful  investigation  of  the  forms  of  the  functions  showed  that  these  numerical 
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infinities  were  removable  by  other  numerical  infinities.  For  example,  an  exponential  in  the 
numerator  and  denominator  of  a  function  would  give  rise  to  two  overflows  when  evaluated 
by  a  computer  but  would  cancel  each  other  and  give  rise  to  a  finite  result  when  analytic 
evaluation  was  properly  applied.  The  construction  of  the  function,  FUNZ,  in  the  TXYZ 
code  was  based  upon  such  evaluations  and  their  numerical  implementations.  With  the 
passage  of  time  and  the  use  of  the  TXYZ  code,  one  place  where  such  cancellation  was 
found  to  be  incomplete  was  in  the  case  of  a  fairly  thick  middle  layer.  For  most  cases 
of  interest  to  the  modeling  of  the  steady-state  thermal  response  of  semiconductor  device 
structures,  these  layers  are  usually  thin.  However,  this  situation  has  now  been  carefully 
investigated  and  the  numerically  stable  limiting  form  has  been  determined  and  used  in  the 
code  in  the  middle-layer  portion  of  the  FUNZ  function.  The  introduction  of  the  nonsingular 
Hmiting  form  is  contained  in  the  TXYZ20  code. 

The  original  TXYZ  program  and  the  TXYZ20  update  have  been  used  for  a  number  of 
investigations.  Some  of  these  include:  the  thermal  evaluation  of  VLSI  packages  using  test 
chips  [8],  the  understanding  of  thermal  resistance  measurements  [9],  the  investigation  of 
the  thermal  interaction  between  electromigation  test  structures  [10],  and  the  benchmarking 
of  other  codes  as,  for  example,  in  the  modeling  of  MMIC  devices  for  the  determination  of 
the  channel  temperatures  during  life  tests  [11].  The  codes  have  been  run  on  a  variety  of 
machines.  They  were  originally  programmed  and  run  on  a  VAX  11/785.  They  now  are 
run  on  PCs  and  Sun  SPARC10  Workstations  [12]. 

PART  2.  ANALYTIC  EVALUATION  OF  LINE  AND  AREA  AVERAGES 

2.1  -  INTRODUCTION 

Electrical  and  optical  measurements  are  often  used  to  determine  the  local  steady-state 
temperature  of  semiconductor  device  structures.  For  example,  the  electrical  resistance  of 
a  segment  of  a  metallization  stripe  may  be  used  to  measure  the  average  temperature  of 
that  segment  of  the  stripe  [13].  Measurements  of  the  temperature  of  a  metal  fine  on  a  thin 
layer  of  silicon  dioxide  deposited  on  a  silicon  substrate  have  been  performed  previously. 
These  have  been  used  in  conjunction  with  the  calculation  of  the  average  temperature  of 
the  metal  fine  in  order  to  determine  the  thermal  conductivity  of  the  silicon  dioxide  layer 
[14]. 

Another  example  is  found  in  an  electrical  technique  for  the  measurement  of  the  peak  junc- 
tion temperature  of  power  transistors  [15].  This  has  recently  been  applied  to  the  use  of 
the  gate  voltage  in  the  measurement  of  the  average  channel  temperature  of  a  power  GaAs 
MESFET  (MEtal- Semiconductor  Field  Effect  Transistor)  [16,17].  The  work  in  reference 
[5]  employs  the  calculated  average  channel  temperature  and  the  measured  average  channel 
temperature  to  determine  a  scale  factor  which  is  used  to  extract  the  peak  channel  tem- 
perature. It  is  the  peak  channel  temperature  which  plays  a  central  role  in  estabhshing 
operating  conditions  for  reliable  device  performance. 

The  temperature  calculations  involved  in  the  above  determinations  of  the  thermal  conduc- 
tivity of  a  Si02  film  [14]  and  the  peak  channel  temperature  of  power  GaAs  MESFETs  [17] 
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make  use  of  the  Kokkas  model  for  the  steady-state  temperature  in  a  rectangular  device 
structure.  In  particular,  the  point  function  temperature,  T(x,t/,  z),  was  calculated  by  us- 
ing the  numerical  implementation  of  the  Kokkas  model  [4]  as  presented  in  the  TXYZ  codes 
[6,7].  However,  as  indicated  in  the  above  discussion,  the  measured  temperature  is  actually 
an  average  over  the  area  over  which  the  electrical  measurement  is  performed.  In  order 
to  make  the  connection  between  the  measured  average  temperature  and  the  calculated 
average  temperature,  references  [14]  and  [17]  use  a  set  of  assumed  representative  points 
at  which  T(x,y,z)  is  calculated.  These  values  are  then  used  to  construct  what  might  be 
called  a  point-by-point  evaluation  of  the  average. 

A  number  of  questions  need  to  be  addressed  in  this  point-by-point  evaluation  of  the  aver- 
age. First,  for  a  given  heat  source  or  active  area  geometry,  how  does  one  choose  a  set  of 
representative  points?  Second,  how  many  points  are  sufficient  to  adequately  represent  the 
effects  of  variations  of  T(x,  y,  z)  over  the  area?  Third,  what  are  the  effects  of  any  residual 
oscillations  in  T(x,y,z)  which  are  inherent  in  the  Fourier  analysis  involved  in  the  model? 

Fortunately,  the  Kokkas  model  and  the  TXYZ  numerical  implementation  provide  for  an 
analytic  evaluation  of  the  averages  over  arbitrary  areas  as  well  as  fine  segments.  This  gives 
the  averages  directly  and  thus  obviates  consideration  of  questions  related  to  the  choice  of 
a  set  of  representative  points  as  well  as  the  number  of  points.  In  addition,  the  calculation 
yields  a  uniformly  convergent  value  of  the  average,  thus  bypassing  any  residual  oscillations 
of  individual  values  of  T(x,  y,  z).  The  averages  take  less  time  to  compute  than  the  point- 
by-point  evaluations  and  provide  a  better  representation  of  the  measured  averages. 


2.2  -  CALCULATION  OF  AVERAGE  TEMPERATURES 
From  eq  (95),  the  point  function  temperature,  T{(x,y,z),  is 


rp  /         x      „  V"  \  -  ±U{n,m)Ti{n,m,z)  cos(nnx/Lx)  cos(miry/Ly) 

Ti{x,y,z)  =  PQ  >    V   77 — 7-7777-  —777  T  •  (102) 

{OnO  +  l)(om0  +  l)LxLyK3 


n=0  m = 0 


As  indicated  above,  the  point -by- point  evaluation  of  the  area  average  temperature  makes 
use  of  a  series  of  point  values  given  by  eq  (102). 

Alternatively,  the  average  temperatures  may  be  calculated  directly.  The  development  ap- 
plies to  both  line  and  area  averages.  Line  averages  provide  a  certain  amount  of  information, 
but  it  is  the  area  average  which  is  related  to  the  experimentally  measured  temperature. 
Consider  an  arbitrary  fine  segment  from  (xj)  to  (xj  +  Ixj)  and  an  arbitrary  rectangular 
area  from  (xj,yj)  to  (xj  +  lxj,yj  +  lyj).  These  are  pictured  in  figure  2  as  well  as  in  the 
top  view  shown  in  figure  3.  It  is  important  to  note  that  these  are  shown  on  the  surface 
(z  —  0)  of  the  structure  for  the  purposes  of  clarity  and  illustration.  The  arbitrary  fine 
segment  and  arbitrary  area  (in  the  x-y  plane)  may  be  located  on  or  inside  the  structure 
(0  <  z  <  -(Lj  +  L2  +  Ls)). 
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Line  for  average 


TOP  VIEW 

Figure  3.  Top  view  of  the  surface  of  the  structure  shown  in  figure  2.  The  line  for  averaging 
and  area  for  averaging  are  separated  from  the  heat  source  for  the  purpose  of  clarity. 
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The  z- dependence  is  explicitly  expressed  in  terms  of  the  averages  which  are  defined  below. 
The  average  of  T{(x,  y,  z)  over  the  arbitrary  line  segment  is  defined  as 


{Ti(y,z))i  =  ~  I  3+  3  Ti(x,y,z)dx,  (103) 
and  the  average  of  T{(x,y,z)  over  the  arbitrary  rectangular  area 

(Ti(z))a  =  -  —  /  /  Tt(x,y,z)dxdy.  (104) 

ix j  ■  lyj  j x.  jy. 

Notice  that  the  x  and  y  dependence  of  Ti(x,y,z)  is  completely  explicit.  This  means  that 
the  integrals  in  the  line  and  area  averages  may  be  analytically  evaluated.  It  is  convenient 
to  define  the  functions,  Ax(n)  and  Ay(m),  by 

,  .      Lx  I  .   (nitix  j+lxj)  \  /mrxj\  1 

A'(n)-^(sm(   l,   )-sm(-ir)[  <105> 

and 

vm,^{s^r^^)-5I„(^)}.  doc, 

By  substituting  eq  (102)  into  eqs  (103)  and  (104)  and  making  use  of  the  definitions  of  the 
functions,  Ax(n)  and  Ay(m),  given  by  eqs  (105)  and  (106),  it  is  straightforward  to  evaluate 
the  fine  and  area  averages. 

Then,  the  average  of  Ti(x,y,z)  over  the  arbitrary  fine  segment  is 


irnt      w       -fo  4U(n,rn)Ti(n,m,z)Ax(n)cos(rrnry/Ly) 

IXj  ^-J  ^  (6n0  +  l)(0m0  +  l)LxLyK3 


and  the  average  of  Ti(x,  y,  z)  over  the  arbitrary  rectangular  area  is 


(107) 


Pq  4JJ(n,m)Ti(n,m,  z)Ax(n)Ay(m) 

{    l[z))a=JX~ly~JZ^Ql^Q        (6nO  +  l)(6mO+l)LxLyK3       '  ^ 
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Several  observations  about  the  line  and  area  average  temperature  defined  by  eqs  (107)  and 
(108)  are  important  to  note.  First,  the  averages  may  be  calculated  by  a  single  evaluation 
of  the  double  sum.  This  means  that  the  averages  may  be  calculated  about  as  quickly  as 
the  evaluation  of  the  point  value  of  the  temperature,  T;(x,  y,  z).  Second,  the  averages  may 
include  any  part  or  all  of  the  x-y  plane  defined  by  the  structure.  If  the  area  is  on  the 
top  surface,  it  also  may  contain  none,  part,  or  all  of  the  heat  source(s).  This  flexibility 
provides  for  the  evaluation  of  area  average  temperatures  over  a  portion  of  the  heat  source 
or  over  a  portion  of  the  surface  away  from  the  heat  source(s).  Third,  the  averages  may  be 
calculated  for  the  surface  (z  =  0)  as  well  as  any  depth  (z  <  0)  inside  the  structure.  Area 
averages  calculated  for  different  areas  and  as  a  function  of  depth  provide  information  on 
the  temperature  spreading  and  heat  transport  inside  the  structure. 

2.3  -   NUMERICAL  IMPLEMENTATION  AND  TXYZ30 

The  basic  TXYZ  program  has  been  used  as  the  framework  in  which  to  construct  the 
programs  to  numerically  evaluate  the  equations  for  the  line  and  area  averages  in  eqs  (107) 
and  (108).  The  resulting  FORTRAN77  code  has  been  named  TXYZ30.  This  code  has  been 
constructed  to  allow  for  the  calculation  of  point  functions,  line  averages,  or  area  averages 
by  means  of  an  integer  variable  switch  which  is  set  in  the  beginning  of  the  input  data  file. 

It  is  important  to  keep  in  mind  that  the  fines  and  areas  over  which  the  averages  are 
computed  are  arbitrary.  They  may  be  on  the  top  surface  (z  —  0)  or  inside  the  structure 
{z  <  0).  In  addition,  the  averages  may  be  easily  computed  not  only  over  part  or  all  of 
the  heat  sources  but  also  at  locations  away  from  the  heat  sources.  Most  applications  of 
these  codes  would  probably  make  use  of  area  averaging.  Line  averaging  may  be  useful  in 
examining  the  fine  average  temperatures.  Typical  evaluations  of  the  averages  take  about 
the  same  CPU  time  as  is  required  to  calculate  a  single  point  value  of  T(x,y,z). 

The  direct  calculation  of  the  area  average  temperature  vastly  simplifies  the  analysis  of 
the  data  when  investigating  the  thermal  conductivity  of  thin  surface  layers.  The  direct 
calculation  of  the  averages  provides  for  a  faster  and  more  accurate  technique  for  the  case  of 
the  determination  of  the  peak  channel  temperature  from  the  ratio  of  the  measured  average 
channel  temperature  and  the  calculated  average  channel  temperature.  This  is  especially 
relevant  when  there  may  be  a  significant  variation  of  T(x,y,z),  both  across  and  along  the 
active  area.  The  point -by- point  evaluation  of  the  average  channel  temperature  must  take 
this  into  account  to  provide  an  adequate  representation  of  the  variation  to  ensure  suitable 
accuracy.  The  use  of  eq  (108)  is  not  encumbered  by  this  requirement  and  provides  the 
averages  directly. 

The  TXYZ30  code  is  contained  in  the  first  appendix.  The  FORTRAN  source  code  is 
also  contained  in  the  HOTPAC  software  package.  In  addition,  there  are  several  sample 
files.  These  are  txyz30io.l,  txyz30io.2  and  txyz30io.3  where  the  io  refers  to  input/output. 
The  file  extension  refers  to  the  example  of  a  point  function  (1),  a  line  average  (2)  and  an 
area  average  (3).  These  files  contain  the  input  data,  the  annotated  input  data  and  the 
corresponding  output  data  for  each  case.  The  user  can  electronically  cut  out  the  input 
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data  for  the  input.dat  file  needed  to  run  TXYZ30. 

2.4  -   TOTAL  AREA  AVERAGES  AND  THE  THERMAL  RESISTANCE 

An  interesting  example  of  the  use  of  eq  (108)  is  provided  by  the  situation  where  the  average 
is  performed  over  the  entire  top  surface  (z  =  0)  area  of  the  structure.  For  the  case  of  a 
single  heat  source  of  dimensions  Ix  by  ly,  it  is  straightforward  to  show  that  the  average 
over  the  total  area  (average  surface  temperature)  is 

(ftP'^Sr.  (io9) 

which  is  just  the  one- dimensional  (thermal  resistance)  result  scaled  by  the  ratio  of  the  heat 
source  area  and  the  total  area. 

2.5  -   APPLICATIONS  OF  AREA  AVERAGES 

The  evaluation  of  the  averages  is  simple  and  direct  and  easily  performed  for  any  number  of 
arbitrary  fine  segments  or  areas.  The  computed  averages  over  part  or  all  of  the  heat  source 
area  should  give  a  more  direct  representation  of  the  measured  values  and  hence  provide  a 
more  accurate  determination  of  the  thermal  properties  as  well  as  peak  temperatures.  In 
addition,  the  averages  are  not  restricted  to  the  areas  where  heat  enters  the  structure  (over 
part  or  all  of  the  heat  source  or  sources).  The  averages  may  also  be  computed  for  various 
depths  inside  the  structure.  Both  nonheat  source  averages  and  depth  averages  can  provide 
information  on  the  effects  of  structure  geometry  and  thermal  properties  on  the  heat  flow. 

Infrared  imaging  techniques  are  also  used  to  experimentally  measure  the  surface  tempera- 
ture of  semiconductor  device  structures.  The  temperature  is  an  average  over  the  spot  size 
of  the  imager.  Application  of  the  calculation  of  the  area  averages  for  this  situation  will 
speed  the  analysis  of  the  average  temperatures  being  measured.  This  should  substantially 
improve  the  understanding  and  application  of  the  techniques. 

Before  turning  to  the  thermal  multilayer  problem,  it  is  important  to  notice  that  the  ap- 
plication of  the  averaging  technique  is  general  and  depends  upon  the  explicit  x-  and  in- 
dependence of  the  temperature  expressed  in  eq  (102).  The  thermal  multilayer  problem  will 
have  exactly  the  same  kind  of  explicit  x-  and  i/- dependence  so  that  the  averaging  can  be 
applied  directly.  This  is  discussed  in  the  next  part  of  the  report. 

PART  3.  RECURSION  RELATION  SOLUTION  OF  MULTILAYER  MODEL 

3.1  -  INTRODUCTION 

It  has  been  two  decades  since  Kokkas'  original  work  on  the  thermal  analysis  of  multiple- 
layer  structures  [3,4].  This  is  based  on  the  solution  of  the  heat  flow  equation  for  a  multilayer 
rectangular  structure.  While  the  interface  boundary  conditions  for  the  continuity  of  the 


31 


temperature  and  the  heat  flow  are  for  a  general  multilayer  structure,  the  equations  have 
been  solved  in  detail  and  evaluated  for  the  case  of  up  to  a  three-layer  structure.  The 
results  of  the  three-layer  steady-state  heat  flow  analysis  are  exact  for  this  situation.  This 
three-layer  solution  has  been  numerically  implemented  in  the  TXYZ  programs  [3,4].  This 
model  and  the  codes  have  found  wide  use  in  a  variety  of  steady-state  problems. 

During  this  same  20-year  time  span,  there  have  been  a  number  of  advances  which  have 
taken  place  in  the  theoretical  analysis  of  the  analogous  electrical  problem  of  two-probe 
spreading  resistance.  The  starting  point  of  the  work  on  the  electrical  problem  is  the 
multilayer  Laplace  equation  analysis  of  Schumann  and  Gardner  [18,19].  Originally,  the 
numerical  analysis  required  a  mainframe  computer  to  do  the  rather  cumbersome  matrix 
algebra  which  develops  quickly  once  the  calculations  get  past  two  or  three  layers.  This 
limitation  coupled  with  difficulties  involved  in  the  evaluation  of  the  accompanying  oscil- 
latory integrals  served  to  severely  limit  the  applicability  of  the  technique.  However,  a 
major  advance  came  with  the  introduction  of  a  recursion  relation  by  Choo  and  coworkers 
[20]  for  the  construction  of  the  solution  of  the  iV-layer  problem  from  the  solution  of  the 
(N  —  l)-layer  problem.  The  recursion  relation  had  previously  been  introduced  in  the  geo- 
physical literature  [21,22].  The  book  by  Koefoed  [21]  provides  a  comprehensive  discussion 
of  Laplace's  equation  and  its  solution,  the  multilayer  approach  applied  to  geoelectric  resis- 
tance measurements,  and  the  derivation  of  the  recursion  relation.  It  should  be  noted  that 
the  mathematical  description  of  spreading  resistance  and  geoelectric  resistance  are  exactly 
the  same  even  though  the  length  scales  are  vastly  different. 

The  introduction  of  the  recursion  relation  eliminated  the  cumbersome  matrix  manipula- 
tions and  allowed  for  the  relatively  straightforward  algebraic  evaluation  of  the  functions 
required  for  the  multilayer  analysis.  The  form  of  the  recursion  relation  has  been  clarified 
by  Berkowitz  and  Lux  [23]  and  has  provided  one  of  the  cornerstones  of  the  routine  calcu- 
lation of  the  spreading  resistance  from  the  resistivity  structure  as  well  as  the  extraction 
of  the  resistivity  structure  from  the  spreading  resistance.  The  remaining  cornerstone  was 
introduced  by  inventive  quadrature  techniques  which  allows  for  the  rapid  evaluation  of  the 
necessary  integrals.  Calculations  which  would  have  been  prohibitive  if  not  impossible  20 
years  ago  can  now  be  performed  routinely  on  PC-  based  systems.  The  model,  the  numerical 
analysis,  and  the  FORTRAN  codes  have  been  summarized  recently  [24]. 

As  the  electrical  problem  and  the  thermal  problem  are  similar  from  the  mathematical 
viewpoint,  the  question  arose  as  to  the  applicability  of  the  recursion  relation  technique, 
with  possible  and  appropriate  modification,  to  the  steady-state  thermal  problem.  If  this 
could  be  accomplished,  then  it  would  provide  for  an  exact  solution  of  the  steady-state 
surface  temperature  for  a  multilayer  structure  with  any  number  of  layers.  The  present 
work  describes  the  modified  recursion  relation  and  its  application  to  the  thermal  problem. 
The  salient  features  of  the  recursion  relation  are  presented  along  with  a  comparison  of  its 
results  with  those  of  the  Kokkas  three-layer  solution.  The  use  of  the  recursion  relation 
for  extension  of  the  calculations  beyond  three  layers  is  discussed,  with  emphasis  on  the 
small  increase  in  calculation  time  for  each  added  layer.  The  computational  efficiency  of 
the  recursion  relation  is  also  discussed.  This  new  technique  should  prove  extremely  useful 
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in  the  understanding  of  the  thermal  effects  of  multilayer  structures. 

As  the  similarity  to  the  spreading  resistance  problem  is  part  of  the  motivation,  a  brief 
discussion  of  the  analysis  involved  in  this  problem  is  presented.  This  is  followed  by  the 
corresponding  analysis  for  the  multilayer  thermal  problem.  The  connection  of  these  two 
problems  is  then  made,  along  with  a  discussion  of  the  need  for  a  modified  recursion  relation 
for  the  thermal  case.  Its  form  is  then  determined  and  the  thermal  recursion  relation  is 
presented.  Finally,  the  implications  of  this  thermal  recursion  relation  are  discussed. 

3.2  -   THE  N-LAYER  ELECTRICAL  PROBLEM 

For  the  case  of  the  spreading  resistance  problem,  the  multilayer  solution  of  the  Laplace 
equation  provides  for  the  calculation  of  the  resistance  between  the  contacts  on  the  top 
surface  of  a  nonuniform  conductivity  (resistivity)  structure.  The  fundamental  assumption 
of  the  multilayer  analysis  is  that  Laplace's  equation  is  satisfied  in  each  "layer"  in  the 
material.  The  problem  is  set  up  in  cylindrical  coordinates  to  emulate  the  symmetry  of  a 
single  circular  contact  on  the  top  surface  of  the  material.  The  potential  is  assumed  to  be 
independent  of  the  angular  variable  in  this  coordinate  system.  The  Laplace  equation  may 
then  be  written  as 

V2 V(r,  z)  =  —  V(r,  z)  +  -  —  V(r,  z)  +  —  V(r,  z)  =  0,  (110) 
or1  r  or  ozL 

where  V(r,  z)  is  the  potential,  r  is  the  radial  coordinate,  and  z  is  the  depth  coordinate. 
This  equation  may  be  solved  by  means  of  separation  of  variables,  with  the  result  that  a 
particular  solution  is 

V(r,z)  =exp(-Az)J0(Ar)+exp(+Az)Jo(Ar),  (111) 

where  Jo(Xr)  is  the  Bessel  function  and  A2  is  the  separation  of  variables  constant.  The 
boundary  condition  on  the  r  part  of  the  solution  is  that  V(r,  z)  approaches  zero  as  r  tends 
to  infinity.  This  is  satisfied  by  the  above  for  all  values  of  A.  Then,  the  general  solution  is 
an  integral  of  the  particular  solution  with  a  weighting  factor  and  is  of  the  form 

V(r,z)  =  J    {(l  +  6(\))exv(-\z)  +  j>(\)exp(+\z)}j0(\r)d\,  (112) 

where  the  weighting  functions,  0(A)  and  ^(A),  are  determined  from  the  z-dependent  bound- 
ary conditions.  The  above  is  the  general  solution  of  Laplace's  equation  in  cylindrical  co- 
ordinates for  a  single  layer.  This  provides  the  framework  in  which  the  multilayer  solution 
may  be  addressed. 
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Figure  4.  Schematic  representation  of  the  geometry  used  in  the  Schumann  and  Gardner 
multilayer  Laplace  equation  analysis.  In  this  figure,  di,  ai,  and  Vi  are  the  thickness, 
electrical  conductivity,  and  potential,  respectively,  in  the  i-ih  layer. 
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The  depth-dependent  portion  of  the  multilayer  geometry  is  presented  in  figure  4.  Each 
layer  is  described  by  a  thickness  and  electrical  conductivity.  Charge  neutrality  is  assumed 
in  each  layer  so  that  the  potential  in  each  layer  satisfies  Laplace's  equation.  The  general 
one-layer  solution  given  by  eq  (112)  provides  a  convenient  and  useful  basis  for  the  required 
layer  solution.  Then,  the  solution  of  Laplace's  equation  in  the  z'-th  layer  in  an  TV-layer 
structure  may  be  written  as 

Vi(r,z)  =  J    {(l+^(A))exp(-Az)  +  Vz(A)exp(+A^)}j0(Ar)dA.  (113) 

The  boundary  conditions  used  to  solve  the  system  of  equations  (i.e.,  determine  {0t(X)} 
and  {V'j(A)},  i  =  1,  ...,N)  are  provided  by  conditions  on  the  top  surface,  the  intermediate 
interfaces,  and  the  bottom  surface.  On  the  top  surface,  current  flow  takes  place  only 
through  the  probe  which  is  modeled  as  a  circular  plate  of  radius  a.  Then  the  top  surface 
boundary  condition  is  expressed  as 

dVN(r,z) 

-<?N  o          =  J{r):  (114) 

oz 

where  is  the  electrical  conductivity  of  the  top  layer,  J(r)  is  the  probe  current  density 
[18,19],  and  the  derivative  is  evaluated  at  z  =  0.  On  the  bottom  surface,  the  potential  is 
assumed  to  be  well  behaved  and,  more  specifically,  is  assumed  to  approach  zero.  This  is 
usually  expressed  as 


lim  Vi(r,z)  =  0.  (115) 


Equations  (114)  and  (115)  provide  2  of  the  2N  conditions  required  to  solve  the  iV-layer 
system  of  equations.  The  remaining  2(iV  —  1)  conditions  are  provided  by  requirements 
at  the  interfaces  between  the  layers  where  the  potentials  and  the  current  densities  are 
assumed  to  be  continuous  [18,19].  These  are  expressed  as 


Vi(r,z)  =  Vi-1(r,z),  (116) 

and 


dV,(r,z)  dV.-^z) 

=  (U7) 

where  the  functions  and  their  derivatives  are  to  be  evaluated  at  the  interfacial  boundaries. 
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For  the  case  of  an  iV-layer  structure,  the  substitution  of  eq  (113)  into  the  boundary  condi- 
tions given  by  eqs  (114  -  117)  gives  rise  to  a  set  of  2iV  equations  in  2N  unknowns  ({#2(A)}, 
{^i(A)},  i  —  l,...,iV).  The  analytic  solution  of  this  system  of  equations  requires  the  use 
of  the  Cramer's  rule  of  linear  algebra.  Clearly,  this  can  become  rather  tedious  especially 
since  the  expansion  coefficients  are  functions  of  the  continuous  variable,  A.  It  is  possible 
to  show  that  the  potential  on  the  top  surface  of  the  iV-layer  structure  may  be  written  as 
[18,19] 

Vn^°)  =  2^J0   A  dK  (U8) 

where  the  kernel  function,  An(\)  =  1  +  2#/v(A),  depends  upon  the  electrical  conductivities 
and  thicknesses  of  the  "layers"  in  the  multilayer  structure  (through  the  solution  of  the 
above  system  of  simultaneous  equations).  Schumann  and  Gardner  were  able  to  work  out 
the  system  of  equations  for  the  cases  up  to  i  =  3.  Cases  beyond  three  layers  presented  a 
major  stumbling  block. 

The  development  of  the  recursion  relation  by  Koefoed  [21]  and  its  use  by  Choo  et  aJ. 
[20]  provided  a  major  breakthrough  in  the  evaluation  of  eq  (118)  and  effectively  removed 
the  numerical  difficulty  associated  with  matrix  inversion.  The  philosophy  and  utility  of  a 
recursion  relation  is  that  the  kernel  function,  A/v(A),  for  an  ./V-layer  structure  can  be  easily 
generated  from  the  kernel  of  an  (N  —  1)  layer  structure  by  means  of  an  algebraic  relation. 
This  algebraic  relation  represents  a  reduction  of  the  matrix  algebra  and  the  subsequent 
manipulations  with  determinants.  The  reader  may  wish  to  refer  to  Koefoed's  book  for  a 
detailed  discussion.  In  particular,  if  the  kernel  is  known  for  the  (N  —  1)  layer  case,  then 
the  kernel  for  the  N  layer  case  is  given  by 

<TN-1u(\)  +  aNAN-1(\) 
AN{A)  =   ■  — —  — ,  (119) 

(TN-l  +  <TN«j{X)AN-i(A) 

where 


=  l-exp(-2A^) 
v  '      1 +exp(-2Ac?N)  v  ' 

and  d,N  is  the  thickness  of  layer  N.  In  practice;,  the  recursion  relation  is  begun  with  the 
one-layer  case  from  which  the  two-layer  kernel  is  generated.  Then  the  three-layer  kernel 
is  generated  from  the  two-layer  kernel.  This  sequence  is  repeated  until  the  iV-layer  kernel 
is  determined.  Notice  how  the  conductivity  and  thickness  of  the  iV-th  layer  enter  in  the 
recursion  relation. 

Berkowitz  and  Lux  [23]  have  shown  that  the  recursion  relation  for  the  kernel  function  could 
be  expressed  as 
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Figure  5.  Schematic  representation  of  the  geometry  used  in  multilayer  analysis  of  the 
steady-state  heat  flow  problem.  In  this  figure,  Li,  K{,  and  Tt  are  the  thickness,  thermal 
conductivity,  and  temperature,  respectively,  in  the  2-th  layer. 


37 


AN(\) 


<jn-i  tanh(AJ/v)  +  CAr^TV-i(A) 
cryv-i  +  &N  An -i  (A)  tanh(  Xd^ ) 


(121) 


The  recursion  relation  is  most  commonly  used  in  this  form.  This  is  in  part  due  to  the 
hyperbolic  tangent  function  which  is  bounded  between  (0,1),  making  the  numerical  im- 
plementation relatively  easy.  This  coupled  with  the  quadrature  techniques  introduced  by 
Berkowitz  and  Lux  has  moved  the  analysis  into  everyday  applicability. 


3.3  -   THE  N-LAYER  THERMAL  PROBLEM 


The  depth-dependent  portion  of  the  multilayer  geometry  for  the  thermal  problem  is  shown 
in  figure  5.  Figures  4  and  5  are  presented  to  reinforce  the  connection  between  the  electrical 
and  thermal  problems  and  the  subsequent  development  of  the  thermal  recursion  relation. 

From  the  development  of  Part  1,  it  is  possible  to  show  that  the  temperature  on  the  top 
surface  of  the  iV-layer  structure  may  be  written  as 


rr  i       n\      p              4U (n,  m)TN(n,  m,  0)  cosjnwx  Lx )  cos(miry/Ly ) 
TN{x,y,0)  =  P0  2^   71 — 7~7vx — .i\r  r  1  ^l22> 

t^On^O  (6n0  +l)(dm0  +l)LxLyKNJ 

The  key  elements  in  eq  (122)  are  the  Fourier  coefficients,  tjv(t),  which  depend  upon  the 
thermal  conductivities  and  thicknesses  of  the  layers  in  the  multilayer  structure  (through 
the  solution  of  the  above  system  of  simultaneous  equations). 

As  discussed  previously,  the  solution  of  eq  (122)  for  cases  beyond  three  layers  presented 
a  major  difficulty.  The  analogous  difficulty  with  the  solution  of  the  multilayer  Laplace 
equation  and  its  removal  by  means  of  the  recursion  relation  provides  the  impetus  for  the 
introduction  of  a  recursive  scheme  in  the  solution  of  multilayer  thermal  equation. 

It  is  important  to  note  that  the  factor  of  7  has  been  explicitly  written  in  the  denominator  of 
eq  (122).  This  has  been  done  to  simplify  the  discussion  of  the  recursion  relation  described 
in  the  next  section.  This  recursion  relation  will  be  used  to  take  into  account  the  multilayer 
thicknesses  and  thermal  conductivities. 


3.4  -   THE  THERMAL  RECURSION  RELATION 


The  strong  resemblances  of  eqs  (22)  and  (114),  eqs  (24)  and  (116),  and  eqs  (25)  and  (117) 
lead  to  the  possibility  of  expressing  the  Fourier  expansion  coefficient  for  the  surface  tem- 
perature of  an  iV-layer  structure  in  the  form  of  eq  (121)  with  the  appropriate  transcription 
from  electrical  conductivities  to  thermal  conductivities,  A  to  7,  An(X)  to  77^(7),  and 
to  Ln-  Then  the  thermal  recursion  relation  would  take  the  form 
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kn-1  tanh^Lw)  +  «jyrjy-i(7) 
kn-1  +  ^NTN-iil)  tanh(jLN) ' 


(123) 


As  the  electrical  and  thermal  bottom-layer  boundary  conditions  are  of  different  form,  it  is 
necessary  to  evaluate  the  corresponding  ^(7).  The  form  of  ^(7)  may  be  obtained  from 
eq  (92)  by  setting  2  =  0  and  then  setting  the  thickness  of  the  two  top  layers  equal  to  zero. 
The  result  is  that 


Equations  (123)  and  (124)  are  the  central  results  of  the  determination  of  the  thermal 
recursion  relation.  They  provide  for  the  determination  of  the  Fourier  coefficients  of  the 
surface  temperature  for  an  TV-layer  structure  from  the  set  of  thicknesses  and  thermal  con- 
ductivities (Ln  Kj-,  i  =  1,  TV)  by  repeated  use  of  the  thermal  recursion  relation.  This  process 
begins  with  the  one-layer  Fourier  coefficient  as  given  by  eq  (124).  This  is  substituted  into 
eq  (123)  to  generate  the  two- layer  Fourier  coefficient.  The  two-layer  Fourier  coefficient  is 
substituted  into  eq  (123)  to  determine  the  three-layer  Fourier  coefficient.  This  process  is 
repeated  until  the  Fourier  ('06151  cient  for  the  desired  number  of  layers  is  generated.  Note 
that  there  is  no  restriction  on  the  number  of  layers.  This  means  that  the  exact  steady-state 
surface  temperature  can  now  be  determined  for  a  multilayer  structure  with  an  arbitrary 
number  of  layers.  It  is  important  to  note  that  this  exact  steady-state  surface  temperature 
satisfies  the  heat  flow  equation  as  well  as  the  necessary  boundary  conditions.  It  does  not 
need  verification  with  existing  layered  solutions. 

The  recursion  relation  technique  should  complement  and  possibly  extend  some  of  the  other 
efforts  [25-29]  aimed  at  problems  beyond  three  layers. 


It  is  instructive  to  compare  and  connect  the  results  of  the  thermal  recursion  relation  with 
those  presented  in  the  literature.  The  principal  one  is  based  upon  the  case  of  a  three- 
layer  structure  where  the  result  of  the  thermal  recursion  relation  is  compared  with  that 
obtained  from  the  Kokkas  model  and  the  TXYZ  code.  Using  eqs  (123)  and  (124),  it  is 
straightforward  to  show  that  the  Fourier  coefficient  for  the  three-layer  structure  is  given 


ti(7)  =  tanh(7£i). 


(124) 


3.5  -   RELATION  TO  PREVIOUS  SOLUTIONS 


by 


Mi) 


K2K3  tanh(7Li)  +  kiK3  tanh(7X2)  +  «i«2  tanh(7Ls)  +  k\  tanh(7Li)  tanh(7L2)  tanh(7L3) 
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Figure  6.  Example  of  the  small  increase  in  computation  time  associated  with  increasing 
the  number  of  layers  in  the  thermal  structure.  The  slope  represents  about  a  5-percent  per 
layer  increase. 
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It  is  also  straightforward  to  show  that  this  is  in  exact  formal  agreement  with  the  Kokkas 
and  TXYZ  expressions  for  z  =  0  for  the  three-layer  case.  This  does  not  serve  as  a  limited 
verification  of  the  recursion  relation  but  rather  shows  that  the  two  are  in  agreement,  as 
expected.  In  addition  to  this  exact  formal  agreement,  numerical  simulations  have  been 
performed  with  the  TXYZ  code  and  a  preliminary  code  based  upon  the  recursion  relation. 
As  expected,  these  are  in  numerical  agreement.  One  interesting  result  is  that  the  recursion 
relation  code  actually  runs  about  15  percent  faster.  This  is  due  to  the  more  compact 
nature  of  the  recursion  relation. 

It  is  interesting  to  note  in  eq  (125)  that  if  kx  =  k2  =  «3,  then  ts(j)  is  independent  of 
the  thermal  conductivities  and  reduces  to  tanh(7(Xi  +  L2  +  L$))  which  is  the  result  for  a 
single  layer  of  thickness  equal  to  L\  +  L2  +  L$.  This  line  of  reasoning  can  be  extended  to 
the  iV-layer  equation.  For  the  case  where  the  thermal  conductivities  of  the  layers  become 
equal,  repeated  use  of  the  addition  formula  for  the  hyperbolic  tangent  shows  that  the  form 
of  the  Fourier  coefficient  becomes  tanh{j(Li  +  ■  — h  Ln))  as  expected  for  a  single  layer  of 
thickness  equal  to  L\  +  •  •  •  -f  Ln- 

3.6  -   NUMERICAL  IMPLEMENTATION  OF  THERMAL  MULTILAYER  CODE,  TML 

The  numerical  implementation  of  the  thermal  recursion  relation  solution  is  contained  in 
the  Thermal  MultiLayer,  TML,  code.  The  FORTRAN  source  code  is  listed  in  the  appendix 
and  is  also  contained  in  the  HOTPAC  software  package. 

For  the  sake  of  illustrating  some  of  the  features  of  the  code,  the  TML  program  has  been  used 
in  calculations  of  the  surface  temperature  for  structures  of  up  to  seven  layers.  The  case  of 
seven  layers  is  used  to  provide  a  convenient  range  and  does  not  represent  a  limitation  of  the 
recursion  relation.  Any  number  of  layers  is  possible  and  can  be  easily  achieved  on  modern 
computing  systems.  These  temperature  calculations  were  performed  for  the  purpose  of 
testing  the  computational  efficiency  of  the  recursive  scheme.  The  results  are  presented 
in  figure  6  and  were  obtained  on  a  VAX  11/785.  They  are  shown  for  the  purpose  of 
illustration  of  the  utility  of  the  recursion  relation.  There  is  only  about  a  5  percent  increase 
in  computation  time  for  each  layer  added  to  the  structure.  This  is  strong  evidence  for  the 
effectiveness  of  the  algebraic  nature  of  the  recursion  relation. 

These  test  cases  were  also  run  for  the  situation  where  all  the  thermal  conductivities  were 
equal.  In  keeping  with  the  above  discussion,  the  results  reduced  to  those  for  a  single  layer 

of  thickness  equal  to  L\  -\  VL^.  This  provided  further  support  of  the  recursion  relation 

and  its  numerical  implementation. 

For  the  purposes  of  acquainting  the  user  with  the  TML  code,  a  number  of  examples  are 
contained  in  the  files.  The  particular  files,  tmlio.l,  tmlio.2,  tmlio.3,  contain  the  three-layer 
problem  contained  in  the  corresponding  TXYZ30  related  files,  txyz30io.l,  txyz30io.2,  and 
txyz30io.3.  The  TML  code,  like  the  TXYZ30  code,  contains  the  point  function,  line  average 
and  area  average  options.  These  are  illustrated  by  the  above. 

For  the  purpose  of  showing  the  use  of  TML  for  real  multilayer  cases,  the  files  tmlio.120, 
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tmlio.220,  and  tmlio.320  are  included  in  the  software  package.  These  illustrate  the  case 
of  the  TML  calculation  for  a  20-layer  case  in  the  point  function,  line  average  and  area 
average  modes. 

3.7  -   GENERALIZED  ONE-DIMENSIONAL  THERMAL  RESISTANCE 

The  special  case  of  uniform  surface  coverage  by  a  single  heat  source  for  a  three-layer 
structure  has  been  presented  in  section  1.8  of  this  report.  There  it  was  shown  that  only 
the  7  — >  0(n  =  0,  m  =  0)  term  needs  to  be  considered.  This  same  line  of  reasoning  applies 
to  the  full  TV-layer  case.  Investigation  of  the  recursion  relation  for  small  values  of  7  shows 
for  the  case  of  uniform  surface  coverage  by  a  single  heat  source  that 

N  T 

2^(0)  =  P0  £ (126) 

which  is  just  the  generalized  one- dimensional  thermal  resistance  result. 

In  section  2.4,  the  average  surface  temperature  (over  the  entire  top  surface  (z  =  0)  area 
of  the  structure)  was  calculated  for  the  case  of  a  single  heat  source  of  dimensions  Ix  by 
ly.  The  result  is  presented  in  eq  (109)  and  is  just  the  one- dimensional  thermal  resistance 
scaled  by  the  ratio  of  the  heat  source  area  and  the  total  area.  Using  the  analysis  used  to 
obtain  eq  (126),  it  is  easy  to  show  that  the  average  surface  temperature  for  the  N-layer 
structure  is 

,™  /  ^  Ix  ■  ly  v^.  L, 

(TN(0))a  =Po-r-rY,-' 
1  y  t=i 

which  is  just  the  generalized  one-dimensional  thermal  resistance  scaled  by  the  ratio  of  the 
heat  source  area  and  the  total  area. 

3.8  -   APPLICATION  TO  BURIED  OXIDE  STRUCTURE 

The  thermal  recursion  relation  should  be  especially  useful  in  the  understanding  of  the 
effects  of  multilayer  structure  thickness  and  thermal  conductivities  on  the  surface  tem- 
perature. The  possibility  of  using  the  surface  temperature  as  a  probe  of  the  thermal 
conductivities  of  the  structure  is  particularly  intriguing.  An  example  of  the  use  of  this 
technique  is  found  in  the  case  of  the  steady-state  surface  temperature  of  a  buried  oxide 
structure.  Recent  trends  in  semiconductor  fabrication  make  use  of  the  buried  oxide  as  a 
way  of  electrically  isolating  the  top  active  region.  Probably  the  most  popular  version  of 
this  is  found  in  SIMOX  (Separation  by  IMplanted  OXygen)  structures  where  the  buried 
oxide  is  formed  by  high-dose  implantation  of  oxygen  followed  by  high-temperature  anneal- 
ing [30].  The  depth  of  the  buried  oxide  is  controlled  by  the  energy  of  the  incident  oxygen 
ion  beam. 
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Figure  7.  Example  of  the  application  of  the  thermal  recursion  relation  to  the  calculation 
of  the  surface  temperature  of  an  SOI  (Silicon-On-Insulator)  structure.  The  insert  (upper 
left)  depicts  the  structure  where  the  layer  thicknesses  and  heater  dimensions  are  not  to 
scale.  The  buried  SiC>2  is  0.00001  cm  thick.  Curves  A,  B,  and  C  are  for  the  cases  where  the 
surface  layers  of  Si  are  0.00001  cm,  0.00002  cm,  and  0.00003  cm  thick,  respectively.  Notice 
how  the  temperatures  do  not  follow  a  one-dimensional  interpretation  in  keeping  with  the 
two-dimensional  heat  flow  in  the  more  thermally  conductive  surface  silicon  layer. 
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When  processed  correctly,  the  implantation  damage  should  be  removed  and  a  buried  layer 
of  SiO-2  should  be  formed.  Electrical  isolation  is  certainly  achieved  but  the  thermal  con- 
ductivity of  the  buried  oxide  may  impose  a  price  of  thermal  isolation  which  may  give  rise 
to  higher- than- desired  operating  temperatures  in  the  device.  These  kinds  of  structures  are 
prime  candidates  for  the  investigation  of  self-heating  effects  in  multiple-layered  structures 
[31,32]. 

As  an  example  of  the  recursion  relation  technique,  preliminary  calculations  were  performed 
on  a  series  of  buried  oxide  structures.  These  were  done  to  ascertain  the  possibility  of  using 
surface  temperature  measurements  as  a  probe  of  the  thermal  conductivity  of  the  buried 
oxide.  The  basic  structure  was  modeled  as  a  chip  (rectangular  structure)  with  lateral  di- 
mensions of  1.0275  cm  by  0.7737  cm.  The  vertical  structure  consisted  of  a  surface  layer  of 
Si  over  a  buried  layer  of  SiC>2  over  a  substrate  layer  of  Si.  The  buried  oxide  was  0.00001  cm 
thick  and  the  substrate  Si  layer  was  0.03189  cm  thick.  Three  thicknesses  of  the  surface  Si 
layer  were  used  in  the  calculations.  These  included  0.00001  cm,  0.00002  cm,  and  0.00003 
cm,  and  are  typical  of  SIMOX  structures.  The  thermal  conductivities  of  the  surface  and 
substrate  Si  layers  were  taken  as  the  room  temperature  value  of  1.55  W/cm°K  and  the 
thermal  conductivities  of  the  buried  oxide  were  taken  in  the  range  from  0.001  to  0.015 
W/cm°K  .  A  small  resistive  heat  source  on  the  top  surface  was  modeled  as  a  rectangular 
element  of  lateral  dimension  of  0.001005  cm  by  0.0825  cm.  The  temperature  (normalized 
by  the  power  density,  T  /  Pq)  was  calculated  in  the  area  of  the  heat  source  for  the  ranges  of 
surface  Si  thickness  and  buried  oxide  thermal  conductivities.  This  was  done  to  ascertain: 
1)  the  effects  of  the  surface  silicon  thickness  and  2)  the  sensitivity  of  the  calculated  temper- 
ature on  the  buried  oxide  thermal  conductivity.  The  results  are  presented  in  figure  7  where 
the  calculated  surface  temperatures  (normalized  by  Pq)  in  the  area  of  the  small  surface 
heat  source  are  plotted  against  the  inverse  thermal  conductivity  of  the  buried  oxide  for 
the  three  surface  silicon  thicknesses.  These  results  indicate  the  important  dependence  on 
the  thickness  of  the  surface  silicon  layer.  In  addition,  the  calculated  temperature  decreases 
with  increasing  thickness  of  the  surface  silicon  layer.  If  the  heat  flow  were  one-dimensional 
and  described  by  eq  (126),  then  the  calculated  temperature  should  increase  with  increasing 
thickness  of  the  surface  silicon  layer.  The  fact  that  the  trend  is  in  the  opposite  direction 
points  to  the  two-dimensional  heat  flow  caused  by  the  more  thermally  conductive  surface 
silicon  layer.  This  necessitates  careful  interpretation  of  the  temperature  data  on  these 
structures. 

3.9  -  SUMMARY  AND  CONCLUSIONS 

Using  the  strong  mathematical  similarity  of  the  multilayer  Laplace  equation  of  spreading 
resistance  and  the  multilayer  steady-state  heat  flow  equation,  a  recursion  relation  has  been 
developed  for  the  multilayer  thermal  problem.  There  are  no  restrictions  on  the  number  of 
layers  which  can  be  used  in  the  calculations.  This  makes  the  results  truly  multilayer.  This 
recursion  relation  has  been  shown  to  produce  the  Kokkas-TXYZ  results  for  the  three-layer 
case.  It  also  gives  rise  to  a  generalized  one- dimensional  thermal  resistance  result  for  the 
case  of  uniform  surface  coverage.  A  preliminary  program  based  on  this  recursive  scheme 
has  been  shown  to  provide  very  good  computational  speed.  As  an  example,  it  has  been 
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applied  to  investigate  the  possible  determination  of  the  buried  oxide  thermal  conductivity 
in  SIMOX-type  structures.  This  recursion  relation  technique  is  simple,  elegant,  and  pow- 
erful and  should  be  extremely  useful  in  the  investigation  and  understanding  of  the  effects 
of  layer  thicknesses  and  thermal  conductivities  on  the  steady-state  surface  temperatures 
of  multilayer  structures  [33,  34]. 

AVAILABILITY  OF  HOTPAC  SOFTWARE  PACKAGE 

The  source  codes  for  the  HOTPAC  software  package  are  written  in  FORTRAN77  and  have 
been  run  on  a  VAX  11/785  minicomputer  system  as  well  as  a  Sun  SPARC10  workstation. 
These  codes  should  run  on  other  systems  with  little,  if  any,  need  for  modification. 

There  are  two  programs  contained  in  the  HOTPAC  software  package.  The  FORTRAN77 
source  code  and  sample  input  and  output  data  files  are  available  in  ASCII  format  on  DOS- 
formatted  floppy  disks.  This  package  is  self-contained  and  is  straightforward  to  run  once 
the  FORTRAN  is  compiled  and  linked  by  the  user-supplied  software.  The  sample  input 
and  output  data  files  are  included  so  that  the  user  can  check  the  programs  for  proper 
operation  as  well  as  to  become  acquainted  with  the  setup  and  use  of  the  codes. 

For  more  information  or  to  receive  a  copy  of  the  code  and  the  report,  please  contact: 

John  Albers 

NIST,  Bldg.  225,  Room  B-310 
Gaithersburg,  MD  USA  20899 
Telephone:  1-301-975-2075 
FAX  Number:  1-301-948-4081 

Netmail  (Internet)  address:  Albers@sed.eeel.nist.gov 
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TXYZ30  I/O  FILE  LISTING  -  txyz30io.l 


1 

10  10 
1  .1 

.5  1.00 
.4  .5 
50  50 
1 

11  0.0  1. 

1 

1  5.  0.0 
1 

1  0.0  0.0 

1  1.0 

1  4.0     2.0     4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
txyz30  point  function  calculation.     This  input  is  annotated  below. 


INPUT 

DESCRIPTION  OF  INPUT 

1 

itype   (=1  for  point  function  calculation) 

10 

10 

x  and  y  dimensions  of  rectangular  structure 

1 

0.1 

thickness  and  thermal  conductivity  of  top  layer  (3) 

.5 

1.00 

thickness  and  thermal  conductivity  of  middle  layer 

.4 

.5 

thickness  and  thermal  conductivity  of  bottom  layer 

50 

50 

upper  summation  limits  for  n  and  m  summations 

1 

iedgex  (=1,  then  read  in  three  values  on  next  line) 

11 

0.0  1. 

number  of  values,   first  point,  increment 

1 

iedgey  (=1,  then  read  in  three  values  on  next  line) 

1 

5.  0.0 

number  of  values,   first  point,  increment 

1 

iedgez   (=1,   then  read  in  three  values  on  next  line) 

1 

0.0  0.0 

number  of  values,   first  point,  increment 

1 

1.0 

number  of  sources  and  power  density 

1 

4.0     2.0     4.0  2.0 

weight,   x,   length  along  x,   y,   length  along  y 

Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  txyz30  using  the  above  input. 


STEADY-STATE  THERMAL  CALCULATION  KOKKAS  ANALYSIS 
POINT  FUNCTION  EVALUATION  OF  T (X, Y, Z) 


X 

Y 

z 

T  (X,  Y,  Z) 

0 

00000 

5. 

00000 

0 

00000 

0. 

02264 

1 

00000 

5. 

00000 

0 

00000 

0. 

03674 

2 

00000 

5. 

00000 

0 

00000 

0. 

15234 

3 

00000 

5. 

00000 

0 

00000 

0. 

65010 

4 

00000 

5. 

00000 

0 

00000 

4. 

30401 

5 

00000 

5. 

00000 

0 

00000 

7  . 

56167 

6 

00000 

5. 

00000 

0 

00000 

4. 

30401 

7 

00000 

5. 

00000 

0 

00000 

0. 

65010 

8 

00000 

5. 

00000 

0 

00000 

0. 

15234 

9 

00000 

5. 

00000 

0 

00000 

0. 

03674 

10 

00000 

5. 

00000 

0 

00000 

0. 

02264 

10.00 

LY= 

10.00 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L3=         1.00000     L2=        0.50000     Ll=  0.40000 
K3=         0.10000     K2=        1.00000     Kl=  0.50000 
UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
NUMBER  OF  HEAT  SOURCE S=  1 
POWER  DENSITY=  1.000000 

WEIGHTS,    COORDINATES,    LENGTHS,   WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 
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TXYZ30  I/O  FILE  LISTING  -  txyz30io.2 


2 

10  10 
1  .1 

.5  1.00 
.4  .5 

50  50 
1 

11  0.0  10. 

1 

1     5.  0.0 
1 

1     0.0  0.0 
1  1.0 

1  4.0  2.0  4.0  2.0 
1 

4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
txyz30  line  average  calculation.     This  input  is  annotated  below. 


INPUT 

2 

10  10 

1  .1 

.5  1.0 

.4  .5 

50  50 

1 

11  0.0 

1 

1 

1  50. 

0. 

0 

1 

1  0.0 

0. 

0 

1  1.0 

1  4.0 

2. 

0 

1 

4.0  2. 

0 

DESCRIPTION  OF  INPUT 

itype  (=2  for  line  average  calculation) 
x  and  y  dimensions  of  rectangular  structure 
thickness  and  thermal  conductivity  of  top  layer  (3) 
thickness  and  thermal  conductivity  of  middle  layer  (2) 
thickness  and  thermal  conductivity  of  bottom  layer  (1) 
upper  summation  limits  for  n  and  m  summations 
iedgex   (=1,   then  read  in  three  values  on  next  line) 
number  of  values,   first  point,  increment 
iedgey  (=1,  then  read  in  three  values  on  next  line) 
number  of  values,   first  point,  increment 
iedgez   (=1,  then  read  in  three  values  on  next  line) 
number  of  values,   first  point,  increment 
number  of  sources  and  power  density 
4.0     2.0      weight,   x,   length  along  x,   y,   length  along  y 
nline   (=1  for  one  line)  as  itype=2 
x,   length  along  x  for  line  element 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  txyz30  using  the  above  input. 


STEADY-STATE  THERMAL  CALCULATION  KOKKAS  ANALYSIS 
LINE  AVERAGE  EVALUATION  OF  TEMPERATURE 

LINE  #  XLINE  LXLINE  Y  Z  AVE  TEMP 

1  4.00000  2.00000  5.00000  0.00000  6.71515 

LX=       10.00         LY=  10.00 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L3=         1.00000     L2=         0.50000     Ll=  0.40000 
K3=         0.10000     K2=         1.00000     Kl=  0.50000 
UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
NUMBER  OF  HEAT  SOURCES=  1 
POWER  DENSITY=  1.000000 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  LINES=  1 

1  4.00000  2.00000 
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TXYZ30  I/O  FILE  LISTING  -  txyz30io.3 


3 

10  10 
1  .1 

.5  1.00 
.4  .5 
50  50 
1 

11  0.0  10. 
1 

1     5.  0.0 
1 

1     0.0  0.0 
1  1.0 

1  4.0  2.0  4.0  2.0 
1 

4.0     2.0     4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
txyz30  area  average  calculation.     This  input  is  annotated  below. 


INPUT  DESCRIPTION  OF  INPUT 

3  itype   (=3  for  area  average  calculation) 

10  10  x  and  y  dimensions  of  rectangular  structure 

1   .1  thickness  and  thermal  conductivity  of  top  layer  (3) 

.5  1.00  thickness  and  thermal  conductivity  of  middle  layer  (2) 

.4   .5  thickness  and  thermal  conductivity  of  bottom  layer  (1) 

50  50  upper  summation  limits  for  n  and  m  summations 

1  iedgex   (=1,   then  read  in  three  values  on  next  line) 

11  0.0     1.  number  of  values,    first  point,  increment 

1  iedgey   (=1,   then  read  in  three  values  on  next  line) 

1     5.     0.0  number  of  values,   first  point,  increment 

1  iedgez   (=1,   then  read  in  three  values  on  next  line) 

1     0.0     0.0  number  of  values,   first  point,  increment 

1       1.0  number  of  sources  and  power  density 

1     4.0     2.0     4.0     2.0      weight,   x,   length  along  x,   y,   length  along  y 

1  narea   (=1  for  one  area)  as  itype=3 

4.0     2.0     4.0     2.0  x,   length  along  x,   y,   length  along  y  for  area 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  txyz30  using  the  above  input. 


STEADY-STATE  THERMAL  CALCULATION  KOKKAS  ANALYSIS 
AREA  AVERAGE  EVALUATION  OF  TEMPERATURE 

AREA  #  XAREA  LXAREA  YAREA  LYAREA  Z  AVE  TEMP 

1  4.00000  2.00000  4.00000  2.00000         0.00000  5.97551 

LX=       10.00         LY=  10.00 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L3=         1.00000     L2=         0.50000     Ll=  0.40000 
K3=         0.10000     K2=         1.00000     Kl=  0.50000 
UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
NUMBER  OF  HEAT  SOURCE S=  1 
POWER  DENSITY=  1.000000 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  AREAS=  1 

1  4.00000  2.00000  4.00000  2.00000 
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TML  I/O  FILE  LISTING  -  tmlio.l 


1 

10  10 
3 

1  .1 
.5  1.00 
.4  .5 
50  50 
1 

11  0.0  1. 

1 

1  5.  0.0 

1  1.0 

1  4.0     2.0     4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  point  function  calculation.     This  input  is  annotated  below. 

INPUT  DESCRIPTION  OF  INPUT 

1  itype   (=1  for  point  function  calculation) 

10  10  x  and  y  dimensions  of  rectangular  structure 
3  number  of  layers  in  the  structure 

1   .1  thickness  and  thermal  conductivity  of  layer  3 

.5  1.00  thickness  and  thermal  conductivity  of  layer  2 

.4   .5  thickness  and  thermal  conductivity  of  layer  1 

50  50  upper  summation  limits  for  n  and  m  summations 

1  iedgex  (=1,   then  read  in  three  values  on  next  line) 

11  0.0     1.  number  of  values,   first  point,  increment 

1  iedgey  (=1,  then  read  in  three  values  on  next  line) 

1     5.     0.0  number  of  values,   first  point,  increment 

1      1.0  number  of  sources  and  power  density 

1     4.0     2.0  4.0     2.0      weight,   x,   length  along  x,   y,   length  along  y 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input. 


STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 
USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 
POINT  FUNCTION  EVALUATION  OF  SURFACE  T(X,Y) 
X  Y  T(X,Y) 

0.00000  5.00000  0.02264 

1.00000  5.00000  0.03674 

2.00000  5.00000  0.15234 

3.00000  5.00000  0.65010 

4.00000  5.00000  4.30401 

5.00000  5.00000  7.56167 

6.00000  5.00000  4.30400 

7.00000  5.00000  0.65010 

8.00000  5.00000  0.15234 

9.00000  5.00000  0.03674 

10.00000  5.00000  0.02264 

LX=       10.00         LY=  10.00 
NUMBER  OF  LAYERS   IN  STRUCTURE=  3 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L  3=         1.00000     K  3=  0.10000 
L  2=         0.50000     K  2=  1.00000 
L  1=         0.40000     K  1=  0.50000 

UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
POWER  DENSITY=  1.000000 
NUMBER  OF  HEAT  SOURCE S=  1 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 
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TML  I/O  FILE  LISTING  -  tmlio.2 

2 

10  10 
3 

1  .1 

.5  1.00 

.4  .5 

50  50 
1 

11  0.0  1. 

1 

1     5.  0.0 
1  1.0 

1     4.0     2.0     4.0  2.0 
1 

4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  line  average  calculation.     This  input  is  annotated  below. 


INPUT  DESCRIPTION  OF  INPUT 

2  itype   (=2  for  line  average  calculation) 

10  10  x  and  y  dimensions  of  rectangular  structure 

3  number  of  layers  in  the  structure 

1   .1  thickness  and  thermal  conductivity  of  layer  3 

.5  1.00  thickness  and  thermal  conductivity  of  layer  2 

.4   .5  thickness  and  thermal  conductivity  of  layer  1 

50  50  upper  summation  limits  for  n  and  m  summations 

1  iedgex   (=1,   then  read  in  three  values  on  next  line) 

11  0.0     1.  number  of  values,   first  point,  increment 

1  iedgey   (=1,  then  read  in  three  values  on  next  line) 

1     5.     0.0  number  of  values,   first  point,  increment 

1      1.0  number  of  sources  and  power  density 

1     4.0    2.0     4.0    2.0      weight,   x,   length  along  x,  y,   length  along  y 

1  nline   (=1  for  one  line)  as  itype=2 

4.0     2.0  x,    length  along  x  for  line  element 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input. 


STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 
USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 
LINE  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE 
LINE   #  XLINE  LXLINE  Y  AVE  TEMP 

1  4.00000  2.00000  5.00000  6.71515 

LX=       10.00         LY=  10.00 
NUMBER  OF  LAYERS   IN  STRUCTURE=  3 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L  3=         1.00000     K  3=  0.10000 
L  2=         0.50000     K  2=  1.00000 
L  1=         0.40000     K  1=  0.50000 

UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
POWER  DENSITY=  1.000000 
NUMBER  OF  HEAT  SOURCES=  1 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  LINES=  1 

1  4.00000  2.00000 
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3 

10  10 

3 

1  .1 
.5  1.00 
.4  .5 
50  50 
1 

11  0.0  1. 

1 

1     5.  0.0 
1  1.0 

1  4.0  2.0  4.0  2.0 
1 

4.0     2.0     4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  area  average  calculation.     This  input  is  annotated  below. 


INPUT  DESCRIPTION  OF  INPUT 

3  itype   (=3  for  area  average  calculation) 

10  10  x  and  y  dimensions  of  rectangular  structure 
3  number  of  layers  in  structure 

1   . 1  thickness  and  thermal  conductivity  of  layer  3 

.5  1.00  thickness  and  thermal  conductivity  of  layer  2 

.4   .5  thickness  and  thermal  conductivity  of  layer  1 

50  50  upper  summation  limits  for  n  and  m  summations 

1  iedgex   (=1,   then  read  in  three  values  on  next  line) 

11  0.0     1.  number  of  values,   first  point,  increment 

1  iedgey   (=1,   then  read  in  three  values  on  next  line) 

1     5.     0.0  number  of  values,   first  point,  increment 

1      1.0  number  of  sources  and  power  density 

1     4.0     2.0     4.0     2.0       weight,   x,   length  along  x,   y,   length  along  y 

1  narea  (=1  for  one  area)  as  itype=3 

4.0     2.0     4.0     2.0  x,   length  along  x,   y,   length  along  y  for  area 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input. 

STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 

USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 

AREA  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE 

AREA  #  XAREA  LXAREA  YAREA  LYAREA       AVE  TEMP 

1  4.00000  2.00000  4.00000  2.00000  5.97551 

LX=       10.00         LY=  10.00 
NUMBER  OF  LAYERS   IN  STRUCTURE=  3 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 
L  3=         1.00000     K  3=  0.10000 
L  2=         0.50000     K  2=  1.00000 
L  1=         0.40000     K  1=  0.50000 

UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
POWER  DENSITY=  1.000000 
NUMBER  OF  HEAT  SOURCE S=  1 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.000.00  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  AREAS=  1 

1  4.00000  2.00000  4.00000  2.00000 
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1 

10  10 
20 


i  a  a  a  a  a 
.  lUUUUU 

1  A  A  A  A  A 
. lUUUUU 

i  a  a  a  a  a 
.  lUUUUU 

T  T  O  O  A  O 

i  a  a  a  a  a 
. lUUUUU 

1  O  C  O  Q  ^ 

i  a  a  a  a  a 
.  1UUU  UU 

,  141z  D4 

i  a  a  a  a  a 

.100000 

T  C  O  A  O  A 

.100000 

.177828 

.100000 

. 199526 

.100000 

.  223872 

.100000 

.251189 

.100000 

OAT  O  O  O 

.281838 

-1  r\  A  A  A  A 

. 100000 

O  1  £  o  o  o 

.  316228 

.100000 

.354813 

i  r\  r\  r\  r\  r\ 
. 100000 

. i ybio  / 

.100000 

.446683 

.100000 

.501187 

.100000 

.562341 

i  r\n Ann 
B lUUUUU 

. Do  u  y O  / 

1  A  A  A  A  A 

-  lUUUUU 

.  /  U  /  y  1 3 

1   /"\  A  A  A  A 

.100000 

n  n  a  o  o  o 

.100000 

.891250 

.0  1. 

0.0 

.0 

0     2.0  4 

0  2.0 

The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  point  function  calculation  for  a  20  layer  structure. 
This  input  is  annotated  below. 


INPUT  DESCRIPTION  OF  INPUT 

1  itype   (=1  for  point  function  calcu 

10  10  x  and  y  dimensions  of  rectangular 

20  number  of  layers  in  the  structure 

.100000  .100000  thickness  and  thermal  conductivity 

.100000  .112202  thickness  and  thermal  conductivity 

.100000  .125893  thickness  and  thermal  conductivity 

.100000  .141254  thickness  and  thermal  conductivity 

.100000  .158489  thickness  and  thermal  conductivity 

.100000  .177828  thickness  and  thermal  conductivity 

.100000  .19952  6  thickness  and  thermal  conductivity 

.100000  .223872  thickness  and  thermal  conductivity 

.100000  .251189  thickness  and  thermal  conductivity 

.100000  .281838  thickness  and  thermal  conductivity 

.100000  .316228  thickness  and  thermal 

.100000  .354813  thickness  and  thermal 

.100000  .398107  thickness  and  thermal  conductivity 

.100000  .446683  thickness  and  thermal  conductivity 

.100000  .501187  thickness  and  thermal  conductivity 

.100000  .562341  thickness  and  thermal  conductivity 

.100000  .630957  thickness  and  thermal  conductivity 

.100000  .7  07  945  thickness  and  thermal  conductivity 

.100000  .7  94328  thickness  and  thermal  conductivity 

.100000  .891250  thickness  and  thermal  conductivity 


lation) 
structure 


conduct  ivit y 
conductivity 


of 

layer 

20 

of 

layer 

19 

of 

layer 

18 

of 

layer 

17 

of 

layer 

16 

of 

layer 

15 

of 

layer 

14 

of 

layer 

13 

of 

layer 

12 

of 

layer 

11 

of 

layer 

10 

of 

layer 

9 

of 

layer 

8 

of 

layer 

7 

of 

layer 

6 

of 

layer 

5 

of 

layer 

4 

of 

layer 

3 

of 

layer 

2 

of 

layer 

1 
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50  50  upper  summation  limits  for  n  and  m  summations 

I  iedgex   (=1,  then  read  in  three  values  on  next  line) 

II  0.0     1.  number  of  values,    first  point,  increment 

1  iedgey   (=1,  then  read  in  three  values  on  next  line) 

1  5.     0.0  number  of  values,   first  point,  increment 

1  1.0  number  of  sources  and  power  density 

1  4.0     2.0     4.0     2.0  weight,   x,   length  along  x,   y,   length  along  y 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input  for  a  20  layer  structure. 


STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 

USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 

POINT  FUNCTION  EVALUATION  OF  SURFACE  T (X, Y) 


X 

Y 

I 

(X,Y) 

0. 

,00000 

5 

.00000 

0. 

01703 

1. 

,00000 

5 

.00000 

0. 

02522 

2. 

,00000 

5 

.00000 

0. 

10517 

3. 

,00000 

5 

.00000 

0. 

42075 

4. 

.00000 

5 

.00000 

3. 

24066 

5. 

.00000 

5 

.00000 

5. 

81674 

6. 

.  00000 

5 

. 00000 

3. 

24066 

7. 

.00000 

5 

.00000 

0  . 

42075 

3, 

.00000 

5 

. 00000 

0. 

10517 

9. 

.00000 

5 

.00000 

0. 

02522 

10, 

.00000 

5 

.00000 

0. 

01703 

10.00 

LY= 

10.00 

NUMBER  OF  LAYERS   IN  STRUCTURE=20 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 


L20= 

0 

.10000 

K20= 

0. 

,10000 

L19= 

0 

. 10000 

K19= 

0. 

,11220 

L18= 

0 

.10000 

K18= 

0. 

,12589 

L17  = 

0 

.10000 

K17  = 

0. 

,14125 

LI  6= 

0 

.10000 

K16= 

0. 

,15849 

L15= 

0 

. 10000 

K15= 

0  . 

.  17783 

L14= 

0 

. 10000 

K14= 

0  , 

.19953 

L13= 

0 

. 10000 

K13= 

0. 

.22387 

L12= 

0 

. 10000 

K12= 

0. 

.25119 

Lll  = 

0 

. 10000 

Kll  = 

0. 

.28184 

L10  = 

0 

.  10000 

K10= 

0  . 

.31623 

L  9= 

0 

.10000 

K  9= 

0. 

.35481 

L  8= 

0 

.10000 

K  8= 

0. 

.39811 

L  7  = 

0 

. 10000 

K  7  = 

0. 

.44668 

L  6= 

0 

.10000 

K  6= 

0. 

.50119 

L  5= 

0 

. 10000 

K  5= 

0. 

.56234 

L  4= 

0 

. 10000 

K  4= 

0 

. 63096 

L  3= 

0 

.10000 

K  3= 

0 

.70794 

L  2= 

0 

.10000 

K  2= 

0 

.79433 

L  1= 

0 

.10000 

K   1  = 

0 

.89125 

UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
POWER  DENSITY=  1.000000 
NUMBER  OF  HEAT  SOURCE S=  1 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 
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2 

10  10 

20 


100000 

. 100000 

100000 

. 112202 

100000 

■  _i_  \y  \j  w  \j  kj 

.  125893 

. 100000 

.  141254 

. 100000 

.  158489 

. 100000 

.  177828 

. 100000 

.199526 

.  100000 

.223872 

. 100000 

.251189 

.  100000 

.281838 

.  100000 

.316228 

. 100000 

.  354813 

.100000 

.398107 

100000 

446683 

.100000 

.501187 

. 100000 

.562341 

.100000 

.630957 

. 100000 

.707945 

. 100000 

.794328 

100000 

. 891250 

.0  1. 

0.0 

.0 

0  2.0 

4.0  2.0 

1 

4.0  2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  line  average  calculation  for  a  20  layer  structure. 
This  input  is  annotated  below. 


INPUT 
2 

10  10 
20 


DESCRIPTION  OF  INPUT 

itype  (=2 

for 

the  line  average  calculation) 

x  and  y  dimensions  of 

rectangular  structure 

number  of 

layers  in  the  structure 

.100000 

.100000 

thickness 

and 

thermal 

conductivity 

of 

layer 

20 

.100000 

.112202 

thickness 

and 

thermal 

conductivity 

of 

layer 

19 

.100000 

.125893 

thickness 

and 

thermal 

conductivity 

of 

layer 

18 

.100000 

.141254 

thickness 

and 

thermal 

conductivity 

of 

layer 

17 

.100000 

.158489 

thickness 

and 

thermal 

conductivity 

of 

layer 

16 

.100000 

.177828 

thickness 

and 

thermal 

conductivity 

of 

layer 

15 

.100000 

.199526 

thickness 

and 

thermal 

conductivity 

of 

layer 

14 

.100000 

.223872 

thickness 

and 

thermal 

conductivity 

of 

layer 

13 

.100000 

.251189 

thickness 

and 

thermal 

conductivity 

of 

layer 

12 

.100000 

.281838 

thickness 

and 

thermal 

conductivity 

of 

layer 

11 

.100000 

.316228 

thickness 

and 

thermal 

conductivity 

of 

layer 

10 

.100000 

.354813 

thickness 

and 

thermal 

conductivity 

of 

layer 

9 

.100000 

.398107 

thickness 

and 

thermal 

conductivity 

of 

layer 

8 

.100000 

.446683 

thickness 

and 

thermal 

conductivity 

of 

layer 

7 

.100000 

.501187 

thickness 

and 

thermal 

conduct  ivit y 

of 

layer 

6 

.100000 

.562341 

thickness 

and 

thermal 

conductivity 

of 

layer 

5 

.100000 

. 630957 

thickness 

and 

thermal 

conductivity 

of 

layer 

4 

.100000 

.707945 

thickness 

and 

thermal 

conductivity 

of 

layer 

3 

.100000 

.794328 

thickness 

and 

thermal 

conductivity 

of 

layer 

2 

.100000 

.891250 

thickness 

and 

thermal 

conductivity 

of 

layer 

1 
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50  50 
1 

11  0.0 

1 

15.  0 
1  1.0 
1  4.0 
1 

4.0  2.0 


0.0 


2.0 


1. 


4.0 


2.0 


upper  summation  limits  for  n  and  m  summations 

iedgex  (=1,  then  read  in  three  values  on  next  line) 

number  of  values,   first  point,  increment 

iedgey  (=1,  then  read  in  three  values  on  next,  line) 

number  of  values,   first  point,  increment 

number  of  sources  and  power  density 

weight,   x,   length  along  x,  y,   length  along  y 

nline  as  itype=2   (=1  for  one  line,  etc) 

x,   length  along  x  for  line  element 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input  for  a  20  layer  structure. 


STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 
USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 
LINE  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE 
LINE  #  XLINE  LXLINE  Y  AVE  TEMP 

1  4.00000  2.00000  5.00000  5.19144 

LX=       10.00         LY=  10.00 
NUMBER  OF  LAYERS  IN  STRUCTURE=2  0 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 


L2  0= 

0. 

.10000 

K20= 

0 

.10000 

L19= 

0. 

.10000 

K19= 

0 

.11220 

L18= 

0. 

.10000 

K18= 

0 

.  12589 

L17  = 

0. 

.10000 

K17  = 

0 

.14125 

LI  6= 

0. 

.10000 

K16= 

0 

.  15849 

L15= 

0. 

.10000 

K15= 

0 

.17783 

1,14= 

0. 

.10000 

K14= 

0 

.19953 

L13= 

0. 

.10000 

K13= 

0 

.22387 

L12= 

0. 

.10000 

K12= 

0 

.25119 

Lll= 

0. 

.10000 

Kll= 

0 

.28184 

L10= 

0. 

.10000 

K10= 

0 

.31623 

L  9= 

0. 

.10000 

K  9= 

0 

.35481 

L  8= 

0. 

.10000 

K  8= 

0 

.39811 

L  7  = 

0. 

.10000 

K  7= 

0 

.44668 

L  6= 

0. 

.10000 

K  6= 

0 

.50119 

L  5= 

0. 

,10000 

K  5= 

0 

.56234 

L  4= 

0. 

.10000 

K  4= 

0 

.63096 

L  3= 

0. 

.10000 

K  3= 

0 

.70794 

L  2= 

0. 

.10000 

K  2= 

0 

.79433 

L  1= 

0. 

.10000 

K  1  = 

0 

.89125 

UPPER  SUMMATION  LIMITS         NUP=       50     MUP=  50 
POWER  DENSITY=  1.000000 
NUMBER  OF  HEAT  SOURCE S=  1 

WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  LINES=  1 

1  4.00000  2.00000 
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3 

10  10 
20 


50  50 
1 

11  0 
1 


t  rt  r\  r\  r\  rt 
lUUUUU 

.lUUUUU 

1  r\  r\  r\  r\  r\ 

100000 

i  i  o  o  n  o 

.  llzzOz 

100000 

.125893 

lUUUUU 

.  141zo4 

i  n  r\  rt  r\  f\ 
1U0U0U 

t  c  o  yi  o  a 

lUUUUU 

.  1  /  /  OZ  O 

i  n  o  n  n  n 
lUUUUU 

i  n  n  n  n  n 
x  u  u  u  u  u 

.  Z  Z  O  O  /  Z 

i nnnnn 

-LUUUUU 

.  Z  3  X  X  O  _7 

i nnnnn 

-LUUUUU 

.  ^  U  1  U JO 

-LUUUUU 

■  J  1  U/iii  o 

i nnnnn 

X  u  u  u  u  u 

•  J J4010 

100000 

.398107 

i nnnnn 
x  u  u  u  u  u 

.  If!  ODOJ 

100000 

.501187 

100000 

.562341 

100000 

.630957 

lUUUUU 

. / u  /  y4o 

lUUUUU 

9  Q  /I  "3  9  O 

i  f\  r\  r*  rs 
lUUUUU 

. ayizou 

0  1. 

0.0 

0 

2.0  4 

0  2.0 

5. 

1. 
4.0 


0     2.0  4.0 


2.0 


The  above  is  the  input  used  in  the  input.dat  file  for  the  example  of  a 
tml  area  average  calculation  for  a  20  layer  structure. 
This  input  is  annotated  below. 


INPUT 
3 

10  10 
20 


. 100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 
.100000 


100000 
,112202 
,125893 
,141254 
,158489 
,177828 
,199526 
,223872 
,251189 
,281838 
.316228 
,354813 
,398107 
,446683 
.501187 
,562341 
.630957 
,707945 
.794328 
,891250 


DESCRIPTION  OF  INPUT 

itype   (=3  for  the  line  average  calculation) 
x  and  y  dimensions  of  rectangular  structure 
number  of  layers  in  the  structure 
thickness  and  thermal  conductivity 


thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity  o 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 
thickness  and  thermal  conductivity 


of 

layer 

20 

of 

layer 

19 

of 

layer 

18 

of 

layer 

17 

of 

layer 

16 

of 

layer 

15 

of 

layer 

14 

of 

layer 

13 

of 

layer 

12 

of 

layer 

11 

of 

layer 

10 

of 

layer 

9 

of 

layer 

8 

of 

layer 

7 

of 

layer 

6 

of 

layer 

5 

of 

layer 

4 

of 

layer 

3 

of 

layer 

2 

of 

layer 

1 
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50  50 
1 

11     0.0  1. 

1 

1     5.  0.0 
1  1.0 

1  4.0  2.0  4.0  2.0 
1 

4.0     2.0     4.0  2.0 


upper  summation  limits  for  n  and  m  summations 

iedgex   (=1,   then  read  in  three  values  on  next  line) 

number  of  values,   first  point,  increment 

iedgey  (=1,  then  read  in  three  values  on  next  line) 

number  of  values,   first  point,  increment 

number  of  sources  and  power  density 

weight,  x,   length  along  x,  y,  length  along  y 

narea  as  itype=3  (=1  for  one  area,  etc) 

x,   length  along  x,   y,   length  along  y  for  area  element 


Below  is  the  output  contained  in  the  output.dat  file  which  is  calculated 
and  written  by  tml  using  the  above  input  for  a  20  layer  structure. 


STEADY-STATE  THERMAL  MULTILAYER  CALCULATION 

USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS 

AREA  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE 

AREA  #  XAREA  LXAREA  YAREA  LYAREA       AVE  TEMP 

1  4.00000  2.00000  4.00000  2.00000  4.64315 

LX=       10.00         LY=  10.00 
NUMBER  OF  LAYERS   IN  STRUCTURE=20 

THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS 


L2  0= 

0 

10000 

K20= 

0 

10000 

LI  9= 

0 

10000 

K19= 

0 

11220 

L18= 

0 

10000 

K18= 

0 

12589 

L17= 

0 

10000 

K17= 

0 

14125 

LI  6= 

0 

10000 

K16= 

0 

15849 

L15= 

0 

10000 

K15= 

0 

17783 

L14= 

0 

10000 

K14= 

0 

19953 

L13= 

0 

10000 

K13= 

0 

22387 

L12= 

0 

10000 

K12= 

0 

25119 

Lll  = 

0 

10000 

Kll= 

0 

28184 

L10= 

0 

10000 

K10= 

0 

31623 

L  9= 

0 

10000 

K  9= 

0 

35481 

L  8= 

0 

10000 

K  8= 

0 

39811 

L  7  = 

0 

10000 

K  7= 

0 

44668 

L  6= 

0 

10000 

K  6= 

0 

.50119 

L  5= 

0 

10000 

K  5= 

0 

56234 

L  4= 

0 

10000 

K  4= 

0 

63096 

L  3= 

0 

10000 

K  3= 

0 

70794 

L  2= 

0 

10000 

K  2= 

0 

79433 

L  1= 

0 

10000 

K  1= 

0 

89125 

UPPER 

SUMMATION 

LIMITS 

NUP= 

POWER 

DENSITY= 

1.000000 

NUMBER  OF 

HEAT  SOURCE S= 

1 

50    MUP=  50 


WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES 

HEAT  SOURCE         WTSOUR  XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  4.00000  4.00000  2.00000  2.00000 

NUMBER  OF  AR£AS=  1 

1  4.00000  2.00000  4.00000  2.00000 
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q* ***************************************************************************** 

C  TXYZ30  -  TXYZ  VERSION  3.0     -  VERSION  DATE  06/27/95 

C 

C  THIS  IS  VERSION  3 . 0  OF  THE  TXYZ  THERMAL  ANALYSIS  CODE.     THIS  VERSION 

C  RETAINS  THE  CHANGES  MADE  IN  VERSION  2.0    (ASSIGNMENT  OF  ARBITRARY  WEIGHTS, 

C  BOTH  +  AND  -,    TO  THE  HEAT  SOURCES  ON  THE  SURFACE  AND  THE  FIX  OF  A  PROBLEM 

C  WITH  THICK  MIDDLE  LAYERS)   AND  INCORPORATES  THE  FOLLOWING  ADDITIONS: 

C  1)   ANALYTIC  AVERAGING  OF  THE  TEMPERATURE  OVER  ARBITRARY  LINE  SEGMENT (S), 

C  2)   ANALYTIC  AVERAGING  OF  THE  TEMPERATURE  OVER  ARBITRARY  AREA ( S ) , 

C  3)    RENUMBERING  OF  LAYERS  TO  BE  IN  KEEPING  WITH  THE  NOTATION  TO  BE  USED 

C  IN  THE  MULTILAYER  THERMAL  PROBLEM    (SEE  TML  CODE) . 

C  THE  LINE  AND  AREA  AVERAGING  ARE  PERFORMED  ANALYTICALLY  FOR  UNIFORM 

C  WEIGHTING  OF  THE  AVERAGE.      THE  AREA  AVERAGE  SHOULD  PROVIDE  A  CLOSER 

C  LINK  WITH  EXPERIMENT  WHERE  THE  MEASUREMENT  DOES  AN  AVERAGING  OVER 

C  SOME  AREA  ELEMENT    (ELECTRICALLY  OR  OPTICALLY) . 
C 

C  THIS  PROGRAM  CALCULATES  THE  STEADY-STATE  TEMPERATURE  DISTRIBUTION 

C  FOR  A  RECTANGULAR  THREE  LAYER  STRUCTURE  WITH  AN  ARBITRARY  NUMBER 

C  OF  RECTANGULAR  HEAT  SOURCES/SINKS  ON  THE  TOP  SURFACE. 

C  THE  CALCULATION  FOLLOWS  FROM  THE  INPUT  OF  THE  THICKNESSES  AND  THERMAL 

C  CONDUCTIVITIES  OF  THE  THREE  LAYERS. 

C  THE  TEMPERATURE    (RELATIVE  TO  THAT  OF  THE  BOTTOM  HEAT  SINK)    MAY  BE 

C  CALCULATED  AS  A  POINT  FUNCTION,    A  LINE  AVERAGE,    OR  AN  AREA  AVERAGE 

C  ANYWHERE  ON  OR  INSIDE  THE  THREE-LAYER  STRUCTURE. 

C  IT  IS   IMPORTANT  TO  EMPHASIZE  THAT  THE  CALCULATION  IS  GENERAL  FOR  THE 

C  THREE  LAYER  STRUCTURE  AND  THE  APPLICATION  TO  SEMICONDUCTOR  STRUCTURES 

C  IS  A  SPECIAL  CASE. 

C  THE  STARTING  EQUATIONS  USED  ARE  GIVEN  IN  EQUATIONS    (13)-(23),   WITH  S=0 

C  (ZERO  FREQUENCY,    STEADY-STATE  CONDITION),    IN  THE  PAPER  BY  KOKKAS  (BELOW). 

C  THE  LINE  AND  AREA  AVERAGES  FOLLOW  FROM  ANALYTICAL  EVALUATION  OF  THESE 

C  EQUATIONS . 
C 

C  REFERENCES:    THE  ORIGINAL  MATHEMATICAL  ANALYSIS  OF  THE  THREE-LAYER 

C  STRUCTURE  WAS  PERFORMED   IN  THE  PAPER  "THERMAL  ANALYSIS 

C  OF  MULTIPLE-LAYERED  STRUCTURES"  BY  ACHILLES  G.  KOKKAS, 

C  IEEE  TRANS.   ELEC .   DEV.   VOL.   ED-21,   NO.    11,    674-681  (1974). 

C  THIS  PAPER  WAS  DRAWN  FROM  HIS  PHD  THESIS:   A.   G.  KOKKAS, 

C  "ANALYSIS  AND  DESIGN  OF  ELECTROTHERMAL  INTEGRATED  CIRCUITS, " 

C  PH.D.    THESIS,    MIT,  1972. 

C 

C  THE  ORIGINAL  FORTRAN  IMPLEMENTATION  OF  THE  STEADY  STATE 

C  KOKKAS  EQUATIONS  IS  CONTAINED  IN  THE  TXYZ  CODE  AND 

C  WAS  PRESENTED   IN  THE  REPORT   "SEMICONDUCTOR  MEASUREMENT 

C  TECHNOLOGY:    TXYZ:   A  PROGRAM  FOR  SEMICONDUCTOR  IC  THERMAL 

C  ANALYSIS"     BY  JOHN  ALBERS,    NBS  SPECIAL  PUBLICATION  400-7  6 

C  (APRIL  1984) . 

C 

C  VERSION  2.0  OF  THE  TXYZ  CODE  WAS  PRESENTED  IN  THE  REPORT 

C  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:   VERSION  2 . 0  OF  THE 

C  TXYZ  THERMAL  ANALYSIS  PROGRAM:   TXYZ2  0"  BY  JOHN  ALBERS, 

C  NIST  SPECIAL  PUBLICATION  400-89    (JUNE  1992) . 
C 

C  VERSION  3.0  OF  THE  TXYZ   CODE  AND  THE  ACCOMPANYING  THERMAL 

C  MULTILAYER  CODES  ARE  DISCUSSED   IN  THIS  REPORT 

C  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:   HOTPAC :  PROGRAMS 

C  FOR  THERMAL  ANALYSIS  INCLUDING  VERSION  3 . 0  OF  THE  TXYZ 

C  PROGRAM,    TXYZ30,   AND  THE  THERMAL  MULTILAYER  PROGRAM,  TML" 

C  BY  JOHN  ALBERS,    NIST  SPECIAL  PUBLICATION  400-96. 

Q*  ********************************************************         ^  ^      ^  *  ^  *  ^  ^  *  ^*  ****  *^  * 
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C  IN  THE  PRESENT  FORM,    THE  PROGRAM  ALLOWS  UP  TO  500  TERMS  TO  BE  INCLUDED 

C  IN  BOTH  THE  N  SUM    (ALONG  X)    AND  THE  M  SUM    (ALONG  Y) .      TO  GO  BEYOND 

C  THIS  NUMBER,    SUBSTITUTE  THE  FOLLOWING  TWO  LINES  WITH  THE  APPROPRIATE 

C  VALUES  OF  NX  AND  MY  FOR  THE  FIRST  TWO  DIMENSION  STATEMENTS  -  MAKE  SURE  TO 

C  REMOVE  THE  COMMENTS  FROM  THE  NEW  LINES  AND  COMMENT  OUT  THE  REPLACED  LINES, 

C  DIMENSION  X( 100) ,    Y(100),    Z(100),  COSYT(MY) 

C  DIMENSION  ARUZER (NX, MY) ,    ARFUNZ (NX,  MY) 

C  ALSO  REPLACE  THE  TESTING  LINE  ABOVE  THE  LINE  LABELLED   110  WITH 
C  IF    (NUP . GT . NX . OR . MUP . GT . MY )    GO  TO  3999 

Q*  ***************************************************************************** 

DIMENSION  X (100) ,    Y(100),    Z(100),  COSYT(500) 
DIMENSION  ARUZER(500,  500) ,    ARFUNZ (500,  500) 
DIMENSION  WTSOUR (50) ,   XSOUR(50),  YSOUR(50) 
DIMENSION  XLINE(30),    XAREA(30),    YAREA (30) 
REAL  LXLINE(30),    LXAREA (30) ,    L YAREA (30) 

REAL  LXSOUR(50),    LYSOUR(50),    K3,    K2,    Kl,    LX,    LY,    L3,    L2,  LI 
COMMON     K3,    K2,    Kl,    LX,    LY,    L3,    L2,  LI 
COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 
C  INPUT  DATA  IS  READ  FROM  10  AND  OUTPUT  GOES  TO  12. 

C  THESE  ARE  WRITTEN  IN  LOWER  CASE.     MANY  OPERATING  SYSTEMS  TAKE 

C  UPPER  AND  LOWER  CASE  AS  EQUIVALENT.      THE  UNIX  OPERATING  SYSTEM 

C  VIEWS   THE  UPPER  AND   LOWER  CASE  NAMED  FILES  AS  DIFFERENT.  UNIX 

C  DEFAULTS   TO  THE  LOWER  CASE  WHICH   IS  USED  HERE. 

open (  unit=10,    file=' input . dat' ,   status  =  'unknown') 
open (  unit=12,   file=' output .dat' ,   status  =  'unknown') 
PI=3 . 14159265 

STEADY-STATE  THERMAL  CALCULATION  KOKKAS  ANALYSIS' ) 
THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS' ) 


2 

FORMAT (IX 

3 

FORMAT (IX 

4 

FORMAT (IX 

5 

FORMAT (IX 

6 

FORMAT (IX 

1   '  MUP=' 

7 

FORMAT (IX 

8 

FORMAT (IX 

9 

FORMAT (IX 

1   ' LXSOUR' 

10 

FORMAT (4X 

11 

FORMAT (IX 

27 

FORMAT (IX 

22 

FORMAT (IX 

31 

FORMAT (IX 

32 

FORMAT ( IX 

33 

FORMAT (IX 

34 

FORMAT (IX 

35 

FORMAT (IX 

40 

FORMAT (IX 

41 

FORMAT (IX 

42 

FORMAT (IX 

43 

FORMAT (6X 

44 

FORMAT (IX 

45 

FORMAT (IX 

l'AVE  TEMP 

46 

FORMAT (IX 

47 

FORMAT (IX 

1    ' L YAREA' 

48 

FORMAT (IX 

11X,F10.5) 

5) 
5) 


'K3=  ',F10.5,'  K2=  ',F10.5,'  Kl=  ',F10 
'L3=  ',F10.5,'  L2=  ' ,F10.5, '  Ll=  ',F10 
'UPPER  SUMMATION  LIMITS  ' , 2X, '  NUP=' , 15, 
15) 

' NUMBER  OF  HEAT  SOURCES=' , 15) 

'WEIGHTS,  COORDINATES,  LENGTHS,  WIDTHS  OF  SOURCES') 
'  HEAT  SOURCE' , 4X, ' WTSOUR' , 8X, ' XSOUR' , 8X, ' YSOUR' , 7X, 
7X,  '  LYSOUR'  ) 

13,  5X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5) 

'POWER  DENSITY=' , Fll . 6) 

'LX=  ',F7.2,3X, '    LY=  ',F7.2) 

F12  .  4,  2X,F12  .  4,  2X,  F12  .  4,  2X,F12  .  4) 

617) 

'NUMBER  OF  LINES=',I3) 
13, 3X,F10.5, 3X,F10.5) 
'NUMBER  OF  AREAS=' , 13) 

13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5) 
'POINT  FUNCTION  EVALUATION  OF  T(X,Y,Z)') 
'LINE  AVERAGE  EVALUATION  OF  TEMPERATURE') 
'AREA  AVERAGE  EVALUATION  OF  TEMPERATURE') 
'  X'  ,  12X,  '  Y'  ,  12X,  '  Z'  ,  9X,  '  T  (X,  Y,  Z)  '  ) 
F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5) 

' LINE  #' , 5X, ' XLINE' , 7X, ' LXLINE' , 8X, ' Y' , 12X, ' Z' , 9X, 
) 

13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3x, flO 
'  AREA  #' , 5X, ' XAREA' , 7X, ' LXAREA' , 8X, ' YAREA'  ,  7X, 
6X, ' Z' , 7x, 'AVE  TEMP' ) 

13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X, FlO. 5, IX, FlO 


5) 


5, 
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51 

FORMAT (IX,  14,  3X,  14 ) 

52 

FORMAT (IX, F10 . 5,  3X,  F10 . 5) 

53 

FORMAT (IX,  11) 

54 

FORMAT (IX, 14, 3X,F10.5, 3X,F10.5) 

55 

FORMAT ( IX,  Fl 0.5) 

56 

FORMAT (IX,  12,  3X,F10.5) 

57 

FORMAT (IX, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X,  F10 

5) 

88 

FORMAT (IX, 'YOUR  UPPER  LIMIT  OF  SUMMATION 

IS  TOO  LARGE.   TRY  AGAIN') 

q* ***************************************************************************** 

 INPUT  SECTION  


=1  FOR  POINT  FUNCTION,    T(X,Y,Z),  EVALUATION 
=2  FOR  LINE  AVERAGE,    <T(Y, Z)>,  EVALUATION 
=3  FOR  AREA  AVERAGE,    <T(Z)>,  EVALUATION 


c 

c 

c 

THE  FOLLOWING 

c 

c 

I TYPE  

c 

c 

c 

c 

LX  

c 

LY  

c 

L3  

c 

K3  

c 

L2  

c 

K2  

c 

LI  

c 

Kl  

c 

NUP  

c 

MUP  

c 

c 

IEDGEX  

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

IEDGEY  

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

THICKNESS  OF  TOP  LAYER 
THERMAL  CONDUCTIVITY  OF  TOP  LAYER 


THICKNESS  OF  BOTTOM  LAYER 
THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 
UPPER  LIMIT  OF  N  SUM,    X  DIRECTION 
UPPER  LIMIT  OF  M  SUM,    Y  DIRECTION 

INDEX  FOR  GENERATING  THE  VALUES  OF  X  TO  BE  USED 
=1   TO  READ  DATA  FOR  FIXED   INCREMENT  X  VALUES 
=2   TO  READ   IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 
IF   IEDGEX=1   THEN  READ  THE  THREE  VARIABLES    (ON  SAME  LINE) 
ILX     THE  NUMBER  OF  X  VALUES  TO  BE  USED 
XI       THE  VALUE  OF  THE  FIRST  POINT  IN  X 
STEPX    (THE   INCREMENT  IN  X) 
IF   IEDGEX=2   THEN  READ  THE  VARIABLE  AND  ARRAY    (ONE  PER  LINE) 
ILX     THE  NUMBER  OF  X  VALUES  TO  BE  USED 
X(I)    THE  ARRAY  OF  X  VALUES    (1=1, ILX) 


=1  TO  READ  DATA  FOR  FIXED  INCREMENT  Y  VALUES 
=2   TO  READ   IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 
IF   IEDGEY=1  THEN  READ  THE  THREE  VARIABLES    (ON  SAME  LINE) 
ILY     THE  NUMBER  OF  Y  VALUES  TO  BE  USED 
Yl       THE  VALUE  OF  THE  FIRST  POINT  IN  Y 
STEPY    (THE   INCREMENT  IN  Y) 
IF   IEDGEY=2  THEN  READ  THE  VARIABLE  AND  ARRAY    (ONE  PER  LINE) 
ILY     THE  NUMBER  OF  Y  VALUES  TO  BE  USED 
Y(I)    THE  ARRAY  OF  Y  VALUES    (1=1, ILY) 


64 


APPENDIX  A  -  TXYZ30  LISTING 


C  IEDGEZ  INDEX  FOR  GENERATING  THE  VALUES  OF   Z  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED   INCREMENT  Z  VALUES 

C  =2  TO  READ  IN  ARRAY  OF  Z  VALUES  OF  NONFIXED  INCREMENT) 

C  IF   IEDGEZ=1   THEN  READ  THE   THREE  VARIABLES    (ON  SAME  LINE) 

C  ILZ     THE  NUMBER  OF   Z  VALUES  TO  BE  USED 

C  Zl       THE  VALUE  OF  THE  FIRST  POINT  IN  Z 

C  STEPZ    (THE   INCREMENT  IN  Z) 

C  IF   IEDGEZ=2  THEN  READ  THE  VARIABLE  AND  ARRAY    (ONE  PER  LINE) 

C  ILZ     THE  NUMBER  OF  Z  VALUES  TO  BE  USED 

C  Z(I)    THE  ARRAY  OF   Z  VALUES    (1=1, ILZ) 

C 

C  NOTE:   ENTER  THE  Z-RELATED  VARIABLES    (Zl,    STEPZ  OR  THE  Z(I)  ARRAY) 

C  AS  ZERO  OR  POSITIVE  QUANTITIES.      THE  PROGRAM  CONVERTS  THE  FINAL 

C  Z(I)    ARRAY  TO  ZERO  OR  NEGATIVE  QUANTITIES  AS  THE  CALCULATION 

C  TAKES  THE  Z  VARIABLE  TO  BE  ZERO  OR  NEGATIVE. 

C 

C  IMPORTANT:    THE  POINT  FUNCTION  EVALUATION  OF  THE  TEMPERATURE  IS 

C  THE  MOST  ELEMENTAL  VERSION.      IN  ORDER  TO  SIMPLIFY  THE 

C  CONSTRUCTION  OF  THE  DATA  FILES  FOR  LINE  AND  AREA 

C  AVERAGES,    THE  PROGRAM  EXPECTS  TO  SEE  THE  ABOVE  IEDGEX, 

C  IEDGEY  AND   IEDGEZ  DATA.      THIS   IS  READ  EVEN  IF   IT   IS  NOT 

C  USED  FOR  THE  AVERAGE  VERSIONS.      HOWEVER,    THE  LINE  AND 

C  AREA  INFORMATION  MAY  THEN  BE  SIMPLY  APPENDED  TO  THE  END 

C  OF  THE  DATA  FILE   IN  ORDER  TO  RUN  THESE  VERSIONS.      SEE  THE 

C  SAMPLE   I/O  FILES  FOR  AN  EXAMPLE  OF  THIS. 

C 

C  NSOUR  NUMBER  OF  HEAT  SOURCES    (UP  TO  50) 

C  P0  POWER  DENSITY 

C 

C  THE  NEXT  NSOUR  LINES  READ  THE  FOLLOWING  INFORMATION  FOR  THE 

C  HEAT  SOURCES    (WITH  ALL  THE  INFORMATION  FOR  EACH  OF  THE  ELEMENTS 

C  ON  A  SINGLE  LINE) 

C 

C  WTSOUR( I) -WEIGHTING  FACTOR  OF   I-TH  SOURCE 

C  (POSITIVE  FOR  SOURCE,   NEGATIVE  FOR  SINK) 

C  XSOUR(I) — X  COORDINATE  OF  ORIGIN  OF  I-TH  SOURCE 

C  LXSOUR (I) -LENGTH  ALONG  X  DIRECTION  OF   I-TH  SOURCE 

C  YSOUR(I) — Y  COORDINATE  OF  ORIGIN  OF   I-TH  SOURCE 

C  LYSOUR (I) -LENGTH  ALONG  Y  DIRECTION  OF   I-TH  SOURCE 

C 

C  IF  ITYPE=1,    THE  POINT  FUNCTION  CALCULATION  CONTINUES  WITH  THE  ABOVE 

C  SET  OF    (X,Y,  Z)  VALUES 

C  IF   ITYPE=2,    THE  LINE  AVERAGE  CALCULATION  READS  THE  FOLLOWING: 

C  NLINE  THE  NUMBER  OF  LINE  SEGMENTS  TO  DO  THE  AVERAGE 

C  THE  NEXT  NLINE  LINES  THEN  READ 

C  XLINE(J),    LXLINE(J) — THE  LOCATION  AND  LENGTH  OF  THE 

C  J-TH  LINE  ELEMENT 

C  IF  ITYPE=3,    THE  AREA  AVERAGE  CALCULATION  READS  THE  FOLLOWING: 

C  NAREA  THE  NUMBER  OF  AREA  SEGMENTS  TO  DO  THE  AVERAGE 

C  THE  NEXT  NAREA  LINES  THEN  READ 

C  XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA (J) -THE  LOCATION 

C  AND  LENGTHS  OF  THE  J-TH  AREA  ELMENT 

***********************************************************  * 
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C  READ   I TYPE    (ANALYSIS  TYPE,    1  FOR  POINT,    2  FOR  LINES,    3  FOR  AREAS) 

READ (10, *) I TYPE 

C  READ  LX  AND  LY    (THE  X  AND  Y  DIMENSIONS  OF  THE  RECTANGULAR  STRUCTURE) 

READ  (10,*)    LX,  LY 

C  READ  L3  AND  K3    (THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  THE  TOP  LAYER) 

READ  (10,*)    L3,  K3 

C  READ   L2  AND  K2    (THICKNESS  AND   THERMAL  CONDUCTIVITY  OF  THE  MIDDLE  LAYER) 

READ  (10,*)    L2,  K2 

C  READ  LI  AND  Kl    (THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  THE  BOTTOM  LAYER) 

READ  (10, * )    LI,  Kl 

C  READ  NUP  AND  MUP    (UPPER  LIMIT  OF  THE   SUMMATION  OVER  THE   INDEX  N  (X-DIR) 

C  UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE   INDEX  M    (Y-DIR) ) 

READ (10, *) NUP, MUP 

C  NUP  AND  MUP  MUST  BE  LESS  THAN  OR  EQUAL  TO  THE  DIMENSIONALITY  OF 

C  ARUZER  AND  ARFUNZ 

C  THE  NEXT  LINE  ALLOWS  VALUES  OF  NUP  AND  MUP  UP  TO  THE  PRESENT  DIMENSION 

C  OF  ARUZER  AND  ARFUNZ.      TO  GO  BEYOND  THIS  VALUE,    THE  USER  SHOULD  EDIT 

C  THE  DIMENSION  STATEMENTS  ACCORDINGLY  AND  THEN  COMMENT  OUT  OR  MODIFY 

C  THE  NEXT  LINE  OF  CODE 

IF    (NUP .GT. 500 .OR. MUP .GT. 500)   GO  TO  3999 

READ (10, *) IEDGEX 

GOTO    (110, 115) IEDGEX 

110  READ    (10, *) ILX, XI, STEPX 
DO   111   1=1, ILX 
X(I)=X1+(I-1) *STEPX 

111  CONTINUE 
GOTO  119 

115  READ  (10,*) ILX 
DO  116  1  =  1,  ILX 
READ (10, *)X (I) 

116  CONTINUE 

119  READ (10,  *) IEDGEY 
GOTO    (120,125)  IEDGEY 

120  READ    (10, *) ILY, Yl, STEPY 
DO  121   1=1, ILY 

Y  (I) =Y1+ (1-1) *STEPY 

121  CONTINUE 
GOTO  12  9 

125  READ  (10,*) ILY 
DO  126  1=1, ILY 
READ (10, *) Y  (I) 

12  6  CONTINUE 

129  READ (10, *) IEDGEZ 
GOTO    (130,135)  IEDGEZ 

130  READ    (10, *) ILZ, Zl, STEPZ 
Zl=-1 . 0*Z1 

STEPZ=-1 . 0*STEPZ 

DO  131   1=1, ILZ 

Z  (I) =Z1+ (1-1) *STEPZ 

131  CONTINUE 
GOTO  139 

135  READ  (10,*)  ILZ 
DO  136  1=1, ILZ 
READ (10, *) Z  (I) 
Z (I)    =  -1 . 0*Z  (I) 

13  6  CONTINUE 
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C  READ   THE  NUMBER  OF  HEAT  SOURCES  AND  THE  POWER  DENSITY 

C  NOTE-POWER  DENSITY  IS  MULTIPLICATIVE  FACTOR  USUALLY  SET  EQUAL  TO  UNITY 

139  READ (10, *)NSOUR,P0 

C  P0   IS  THE  POWER  DENSITY,    ASSUMED  UNIFORM  FOR  ALL  HEATERS 

C  NSOUR  IS  THE  TOTAL  NUMBER  OF  HEATING  ELEMENTS  ON  THE  SURFACE  OF  THE 

C  THE  TOP  LAYER   (UP  TO  50) 

C  THE  NEXT  LOOP  READS   IN  THE  COORDINATES  OF  THE  ORIGIN  OF  THE 

C  HEATING  ELEMENTS  ALONG  WITH  THEIR  LENGTHS  AND  WIDTHS 

C  THE  WEIGHTING  FACTOR  IS  ALSO  ENTERED    (THIS   IS  REAL,  NONINTEGER) 

DO   140   1=1, NSOUR 

READ (10, * ) WTSOUR ( I ) , XSOUR ( I ) ,    LXSOUR ( I )  ,  YSOUR ( I )  ,  LYSOUR(I) 
C  WTSOUR (I)    IS  THE  WEIGHTING  FACTOR  FOR  THE  I-TH  HEATER  ELEMENT 

C  XSOUR (I)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 

C  LXSOUR (I)    IS  THE  LENGTH  OF  THE   I-TH  HEATER  ALONG  THE  X  DIRECTION 

C  YSOUR (I)    IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF   I-TH  HEATER  ELEMENT 

C  LYSOUR(I)    IS  THE  LENGTH  OF  THE  I-TH  HEATER  ALONG  THE  Y  DIRECTION 

140  CONTINUE 

C  THIS   IF. . .THEN. . .ELSE   IF  CONSTRUCTION  IS  USED  TO  READ  IN  THE  DATA  FOR 

C  THE  LINES  OR  AREAS  TO  BE  CONSIDERED  IN  THE  CALCULATION 

C  IF  THE  ANALYSIS   IS  FOR  A  POINT  FUNCTION,    THEN  GO  OUT  OF  THE   IF. THEN. ELSE 

IF    (ITYPE.EQ.l)  THEN 
GOTO   17  0 

C  IF  THE  ANALYSIS   IS  FOR  A  LINE  AVERAGE,    THEN  READ  THE  NUMBER  OF  LINE 

C  SEGMENTS  AND  THEIR  LOCATION 

ELSE  IF    (ITYPE.EQ.2)  THEN 
C  READ  THE  NUMBER  OF  LINE  SEGMENTS 

READ (10, * ) NLINE 

DO  150     J=l, NLINE 
C  READ  THE  ORIGIN  AND  LENGTH  OF  EACH  LINE  SEGMENT 

C  NOTE  THAT  THE  AVERAGE   IS  ALONG  THE  X-DIRECTION  FOR  GIVEN  Y  VALUES 

C  TO  DO  THE  AVERAGE  ALONG  THE  Y-DIRECTION  FOR  GIVEN  X,    SIMPLY  ROTATE 

C  THE  STRUCTURE  BY  90  DEGREES  AND  USE  THE  CORRESPONDING  NEW  X' S  AND  Y' S 

C  XLINE(J)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  J-TH  LINE  ELEMENT 

C  LXLINE(J)    IS  THE  LENGTH  OF  THE  J-TH  LINE  ALONG  THE  X  DIRECTION 

READ (10, * ) XLINE (J) , LXLINE (J) 
150  CONTINUE 

C  IF  THE  ANALYSIS   IS  FOR  A  AREA  AVERAGE,    THEN  READ  THE  NUMBER  OF  AREAS 

C  AND  THEIR  LOCATIONS 

ELSE   IF    (ITYPE.EQ.3)  THEN 
C  READ  THE  NUMBER  OF  AREAS 

READ (10, * ) NAREA 

DO  160  J=l, NAREA 
C  READ  THE  FOLLOWING 

C  XAREA(J)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  J-TH  AREA  ELEMENT 

C  LXAREA(J)    IS  THE  LENGTH  OF  THE  J-TH  AREA  ALONG  THE  X  DIRECTION 

C  YAREA(J)    IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF  J-TH  AREA  ELEMENT 

C  LYAREA (J)    IS  THE   LENGTH  OF  THE   J-TH  AREA  ALONG  THE  Y  DIRECTION 

READ (10, * ) XAREA ( J) , LXAREA (J) , YAREA (J) , LYAREA (J) 
160  CONTINUE 

END  IF 
17  0  CONTINUE 

P04LK  =  4.0   *  P0  /    (  LX*  LY*  K3) 
PILX  =  PI   /  LX 
PILY  =  PI   /  LY 
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Q*  ***************************************************************************** 

c 

C  END  OF  DATA  INPUT  SECTION 

C  BEGIN  CALCULATION  OF  TEMPERATURE 

C  THE  SUBROUTINES  USED   IN  THE  CALCULATION  ARE: 

C  1)    UZERO(N,M)    -  CALCULATES  THE  FOURIER  COSINE  TRANSFORM  OF  THE 

C  FUNCTION,    U(X,Y),    THE  POWER  DENSITY  FUNCTION  FOR  ALL  OF  THE 

C  HEAT  SOURCES. 

C  2)    FUNZ(N,M, Z)    -  CALCULATES  THE  Z-DEPENDENT  PORTION  OF  THE  SUM 

C  REMEMBERING  THAT  THIS   IS  A  FUNCTION  OF  THE  SUMMATION 

C  INDICES    (N, M) . 

Q*  ***************************************************************************** 

C  CALCULATE  THE  FOURIER  COMPONENTS  OF  THE  HEAT  SOURCES,    U (N, M) 

^* ***************************************************************************** 

DO  300  MM=1,MUP 
M  =  MM  -  1 
DO  250  NN=1,NUP 
N  =  NN  -  1 

ARUZER (NN, MM) =UZERO (N, M) 
2  50  CONTINUE 
300  CONTINUE 

q* ***************************************************************************** 
C  END  OF  U(N,M)    CALCULATION  AND  BEGINNING  OF  MAJOR  LOOP  FOR  Z 

Q* ***************************************************************************** 
DO   6000   IZ=1, ILZ 

Q* ***************************************************************************** 
C  CALCULATE  THE   Z  DEPENDENT  PORTION,    I.E.,    FUNZ(N,M, Z) 

(2*  ****************************************************************************  * 

DO  400  MM=1,MUP 
M  =  MM  -  1 
DO  350  NN=1,NUP 
N  =  NN  -  1 

ARFUNZ (NN, MM) =FUNZ (N, M, Z (IZ) ) * ARUZER (NN, MM) 
350  CONTINUE 
400  CONTINUE 

C  THE  FOLLOWING  IF. THEN. ELSE  CONSTRUCTION  OPERATES  ACCORDING  TO  THE 

C  TYPE  OF  ANALYSIS  TO  BE  USED. 

C  FOR  THE  POINT  FUNCTION  ANALYSIS,    THE  3000  LOOP   IS  USED 

C  FOR  THE  LINE  AVERAGE  ANALYSIS,    THE  4000  LOOP   IS  USED 

C  FOR  THE  AREA  AVERAGE  ANALYSIS,    THE  5000  LOOP   IS  USED 

C  THIS  PORTION  DOES  THE  POINT  FUNCTION  CALCULATION 

IF    (ITYPE.EQ.l)  THEN 
C  BEGINNING  OF  POINT  FUNCTION  ANALYSIS 

WRITE  (12, 2) 

WRITE (12, 40) 

WRITE (12, 43) 

DO  3000     IY=1, ILY 

DO  3100  MM=1,MUP 
M  =  MM  -  1 

COSYT (MM) =COS (FLOAT (M) *Y (IY) *PILY) 
310  0  CONTINUE 
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DO  3000   IX=1, ILX 
SUM=0  .  0 

DO  3300  MM=1,MUP 

M  =  MM  -  1 

DO  3200  NN=1,NUP 

N  =  NN  —  1 

NDN=0 

NDM=0 

IF  (N.EQ.O)  NDN=1 
IF    (M.EQ.O)  NDM=1 

TOP  =  ARFUNZ (NN, MM)    *  COS (FLOAT (N) *X (IX) *PILX)    *  COSYT (MM) 

BOTTOM= (NDN+1) * (NDM+1) 

TSUM=TOP /BOTTOM 

SUM=SUM+TSUM 
32  0  0  CONTINUE 
3300  CONTINUE 

TEMP  =  P04LK  *  SUM 

WRITE    (12, 44)X(IX) , Y (IY) , Z (IZ) , TEMP 
3000  CONTINUE 
C  END  OF  POINT  FUNCTION  CALCULATION 

ELSE   IF    (ITYPE.EQ.2)  THEN 
C  THIS  PORTION  DOES  THE  LINE  AVERAGE  CALCULATION 

WRITE (12, 2) 
WRITE (12, 41) 
WRITE  (12,  45) 
DO  4000     IY=1,  ILY 
DO  4010  MM=1,MUP 
M  =  MM  -  1 

COSYT (MM) =COS (FLOAT (M) *Y ( IY) *PILY) 
4  010  CONTINUE 

DO  4000  J=1,NLINE 
SUM=0 . 0 

DO  4200  MM=1,MUP 

M  =  MM  -  1 

DO  4100  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF    (N.EQ.O)  NDN=1 
IF    (M.EQ.O)  NDM=1 
IF(N.EQ.O)   GO  TO  4160 

TERMX=SIN (FLOAT (N) *PI* (XLINE (J) +LXLINE (J) ) /LX) 
1  -SIN (FLOAT  (N) *PI*XLINE (J) /LX) 

TERMX=TERMX*LX/ (FLOAT (N) *PI) 
GO  TO  4165 
4160       TERMX=LXLINE (J) 
4165  CONTINUE 

TOP  =  ARFUNZ (NN, MM)    *  COSYT (MM)    *  TERMX 
EOTTOM= (NDN+1) * (NDM+1) 
TSUM=TOP /BOTTOM 
SUM=SUM+TSUM 
4100  CONTINUE 

TEMP  =  P04LK  *  SUM  /  LXLINE(J) 
42  00  CONTINUE 

WRITE  (12,  46) J,  XLINE (J) , LXLINE (J) , Y (IY) , Z (IZ) , TEMP 
4000  CONTINUE 
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C  END  OF  LINE  AVERAGE  PORTION  OF  THE  CODE 

ELSE   IF    (ITYPE.EQ.3)  THEN 
C  THIS  PART  PERFORMS  THE  AREA  AVERAGE  CALCULATION 

WRITE  (12, 2) 

WRITE  (12, 42) 

WRITE  (12, 47) 

DO  5000  J=l , NAREA 
SUM=0 . 0 

DO  5100  MM=1,MUP 

M  =  MM  —  1 

DO  52  00  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF    (N.EQ.0)  NDN=1 
IF    (M.EQ.0)  NDM=1 
AREA=0 . 0 

IF(N.EQ.O)   GO  TO  5160 

TERMX  =  SIN (FLOAT (N) *PI* (XAREA (J)    +     LXAREA ( J) ) /  LX) 

-  SIN  (FLOAT (N) *PI*XAREA (J) /  LX) 
TERMX=TERMX*   LX/ (FLOAT (N) *PI ) 
GO  TO  5165 

TERMX=  LXAREA (J) 
IF(M.EQ.O)   GO  TO  5164 
TERMY  =  SIN (FLOAT  (M) *PI* (YAREA (J)    +     L YAREA ( J ) ) /  LY) 

-  SIN (FLOAT (M) *P I* YAREA (J) /  LY) 
TERMY= TERMY*   LY/ (FLOAT (M) *PI) 
GO  TO  5166 

TERMY=  L YAREA (J) 
TERMI=TERMX* TERMY 
AREA=TERMI 

TOP  =  ARFUNZ (NN, MM)    *  AREA 
BOTTOM= (NDN+1) * (NDM+1) 
TSUM=TOP /BOTTOM 
SUM=SUM+TSUM 
CONTINUE 

TEMP  =  P04LK  *   SUM  /    (LXAREA ( J) *LYAREA ( J) ) 
CONTINUE 

WRITE (12, 48) J, XAREA (J) , LXAREA (J) , YAREA (J) , L YAREA ( J) , Z (IZ)  ,  TEMP 
5000  CONTINUE 

END  IF 
6000  CONTINUE 

WRITE  (12, 27)    LX,  LY 

WRITE  (12,  3) 

WRITE(12,5)    L3,    L2,  LI 
WRITE  (12,  4)    K3,    K2,  Kl 
WRITE (12, 6)NUP,MUP 
WRITE (12, 7)NSOUR 
WRITE(12,11)  P0 
WRITE  (12, 8) 
WRITE (12, 9) 
DO  3888  I=l,NSOUR 

WRITE(12,  10) I,WTSOUR(I) ,XSOUR(I) , YSOUR(I)  ,    LXSOUR(I),  LYSOUR(I) 
3888  CONTINUE 


1 

5160 
5165 

1 

5164 
5166 

5200 
5100 
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IF    (ITYPE.EQ.l)  THEN 

GOTO  7000 
ELSE   IF    (ITYPE.EQ.2)  THEN 

WRITE  (12, 32)  NLINE 

DO   6100  J=l, NLINE 

WRITE (12, 33) J, XLINE (J) , LXLINE (J) 
6100  CONTINUE 

ELSE  IF    (ITYPE.EQ.3)  THEN 
WRITE (12, 34) NAREA 
DO  6200  J=l, NAREA 

WRITE  (12,  35)  J,  XAREA(J)  ,  LXAREA(J)  ,  YAREA(J)  ,  LYAREA(J) 
6200  CONTINUE 

END  IF 

GO  TO  7  000 
3999  WRITE (12,  88) 
7000  STOP 

END 

Q* ************************************************  ******** 

C  END  OF  THE  MAIN  PROGRAM 
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C  FUNCTION  UZERO  LISTING 

q* ***************************************************************************** 

C  DESCRIPTION  OF  THE  FUNCTION  UZERO (N, M) 

C  THIS  FUNCTION  CALCULATES  THE  DOUBLE  FOURIER  COSINE  TRANSFORM 

C  OF  THE  POWER  DENSITY  FUNCTION,    U (X, Y) .      THIS   IS  THE  TRANSFORM 

C  FOR  ALL  OF  THE  HEAT  SOURCES.      THE  ASSUMPTION  IS  MADE  THAT  THE 

C  POWER  DENSITY  IS  UNIFORM  AND  EQUAL  TO  UNITY  OVER  THE  SURFACE 

C  OF  THE  HEATING  ELEMENTS.      THAT  IS, 

C  U(X, Y)=l      (XSOUR(I)<=X<=XSOUR(I)+LXSOUR(I)  AND 

C  YSOUR(I) <=Y<=YXOUR(I) +LYSOUR(I)  ). 

C  U(X,Y)=0  OTHERWISE. 

C  UNDER  THESE  CONDITIONS,    IT  IS  POSSIBLE  TO  ANALYTICALLY  EVALUATE 

C  THE  DOUBLE   INTEGRAL  FOR  EACH  HEATING  ELEMENT.     AS  THE  HEATING 

C  ELEMENTS  ARE  ASSUMED  TO  BE   INDEPENDENT,    THE  CONTRIBUTION  FROM 

C  EACH  ELEMENT  MAY  BE  ADDED  TO  OBTAIN  THE  U(N,M)    FOR  ALL. 

C 

C  NONUNIFORM  POWER  DENSITIES  MAY  BE  TAKEN  CARE  OF  BY  USING  THE 

C  WEIGHTING  FACTOR  FOR  EACH  ELEMENT    (READ   IN  THE  MAIN  PROGRAM) 
C* ***************************************************************************** 

FUNCTION  UZERO (N, M) 

DIMENSION  WTSOUR( 50),    XSOUR(50),  YSOUR(50) 

REAL  LXSOUR(50),    LYSOUR(50),    K3,    K2,    Kl,    LX,    LY,    L3,    L2,  LI 

COMMON     K3,    K2,    Kl,    LX,    LY,    L3,    L2,  LI 

COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 

PI=3. 14159265 

UZERO=0 . 0 

DO  500  1=1, NSOUR 

IF(N.EQ.O)   GO  TO  100 

TERMX  =  SIN (FLOAT (N) *P I* (XSOUR (I)    +     LXSOUR(I))/  LX) 
1  -  SIN (FLOAT (N) *PI *XSOUR ( I ) /  LX) 

TERMX=TERMX *   LX/ (FLOAT (N) *PI) 
GO  TO  150 
100  TERMX=  LXSOUR (I) 

150  IF(M.EQ.O)    GO  TO  200 

TERMY  =  SIN (FLOAT (M) *PI* (YSOUR (I)   +     LYSOUR(I))/  LY) 
1  -  SIN (FLOAT (M) *PI* YSOUR (I) /  LY) 

TERMY=TERMY*   LY/ (FLOAT (M) *PI) 
GO  TO  250 
200  TERMY=  LYSOUR (I) 

250  TERMI=TERMX* TERMY 

UZERO=UZERO+TERMI *WTSOUR ( I ) 
500  CONTINUE 
RETURN 
END 


72 


APPENDIX  A  -  TXYZ30  LISTING 

q*  ****************************************************************************  * 

C  FUNCTION  FUNZ  LISTING 

C 

C  DESCRIPTION  OF  THE  FUNCTION  FUNZ (N, M, Z) 

C  THIS  FUNCTION  IS  USED  TO  CALCULATE  THE  Z  DEPENDENT  PART  OF  THE 

C  FUNCTION  USING  THE  S=0  VERSIONS  OF  EQUATIONS    (15) -(17)    OF  KOKKAS . 

C  THESE  ARE  USED   IN  CONJUNCTION  WITH  EQUATIONS    (18) -(22)   ALSO  FOR 

C  S=0    (STEADY-STATE  CONDITION) .      THE  SPECIFIC  FORM  OF  FUNZ  IS 

C  DETERMINED  BY  THE  VALUE  OF  Z,    I.E.,    IF  Z  FALLS  IN  THE  TOP,  MIDDLE, 

C  OR  BOTTOM  LAYER  OF  THE  STRUCTURE.      IN  ADDITION,    SPECIAL  CARE  IS 

C  TAKEN  TO  EVALUATE  FUNZ  FOR  THE  CASE  WHERE  GAMMA=0  AS  THESE  ARE 

C  SIMPLE,    BUT  THE  COMPUTER  DOES  NOT  KNOW  HOW  TO  EVALUATE  LIMITS. 

C  IN  ADDITION,    THE  FORM  OF  THE  THREE  SOLUTIONS  HAS  BEEN  CHANGED  TO 

C  GET  RID  OF  THE  ARTIFICIAL  OVERFLOW  PROBLEMS  COMING  FROM  THE 

C  HYPERBOLIC  FUNCTIONS,    COSH  AND  SINH,    FOR  THE  CASES  WHERE  THE 

C  ARGUMENTS  BECOME  LARGE. 
C 

Q* ***************************************************************************** 

FUNCTION  FUNZ(N,M, Z) 

DIMENSION  WTSOUR( 50) ,    XSOUR(50),  YSOUR(50) 

REAL  LXSOUR(50),    LYSOUR(50),   K3,   K2,   Kl,    LX,    LY,    L3,    L2,  LI 

COMMON     K3,    K2,    Kl,    LX,    LY,    L3,    L2,  LI 

COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 

PI=3 . 14159265 

GAMMA=SQRT ( (FLOAT (N) *PI/  LX) **2  +    (FLOAT (M) *PI/  LY)**2) 

VS=GAMMA*  L3 

VC=GAMMA*  L2 

VI=GAMMA*  LI 

VT=GAMMA* (   L3  +  Z) 

VM=GAMMA* (   L3+  L2+Z) 

VB=GAMMA* (   L3+  L2+  Ll+Z) 

BOTl=TANH (VS) *TANH (VI) 

BOTl=BOTl+(   Kl/  K2) *TANH (VS) *TANH (VC) 
BOT2=(   K2/  K3) *TANH  (VI) *TANH (VC)  +  (  Kl/  K3) 
GFUNC=1 . 0/ (BOT1+BOT2) 
AZ=ABS (Z) 

IF    (AZ.GT.    L3)    GO  TO  500 
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**************************************************** 


IF    (GAMMA. EQ. 0.0)    GO  TO  100 
TERMSl=TANH(VI)+(  Kl/  K2)*TANH(VC) 

TERMS2=(  K2/  K3 ) *TANH (VT) * (TANH (VI ) *TANH (VC) + (  Kl/  K2) ) 

TERMS=TERMS1+TERMS2 

IF    (Z.EQ.0.0)   GO  TO  90 

IF    (VS.GT.5.0.AND.VT.GT.5.0)   GO  TO  80 
IF    (VS.LT.5.0)   GO  TO  10 


20  CONTINUE 

IF    (VT.LT.5.0)    GO  TO  30 

C2=0 . 5*EXP (VT) 

GO  TO  40 
30  C2=COSH(VT) 
40  CONTINUE 

FUNZ=GFUNC * TERMS * C 1 * C2 /GAMMA 

RETURN 

80       FUNZ=GFUNC* TERMS *EXP (GAMMA*Z) /GAMMA 
RETURN 

90       FUNZ=GFUNC* TERMS /GAMMA 
RETURN 

100     FUNZ=(  L3+Z)+(  K3/  K2)*  L2+(  K3/  Kl)*  LI 

RETURN 
500     TOTAL=  L3+  L2 

IF    (AZ .GT. TOTAL)   GO  TO  1500 


10 


CI' 
GO 
CI 


2 .0*EXP (-VS) 
TO  2  0 

1.0/COSH (VS) 
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Q* ***************************************************************************** 

C  MIDDLE  LAYER  CALCULATION 

C  THIS  IS  THE  MIDDLE  LAYER  CALCULATION  WHICH  IS  DEFAULTED  TO  IF 

C  Z  FALLS  INTO  THIS  DOMAIN  OF  DEPTHS 

Q*  ****************************************************************************  * 

IF    (GAMMA. EQ. 0.0)    GO  TO  1000 
TERMC=TANH (VI ) + (   Kl/  K2)*TANH(VM) 

IF    (VS.GT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0)    GO  TO  800 

IF    (VS.LT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0)    GO  TO  900 

IF    (VS.LT.5.0)    GO  TO  250 

Cl=2 . 0*EXP (-VS) 

GO  TO  2  60 
250  C1=1.0/COSH(VS) 
2  60  CONTINUE 

IF    (VC.LT.5.0)   GO  TO  270 

C2=2 . 0*EXP (-VC) 

GO  TO  280 
270  C2=1.0/COSH(VC) 
280  CONTINUE 

IF    (VM.LT.5.0)   GO  TO  320 

C3=0.5*EXP (VM) 

GO  TO  330 
320  C3=COSH(VM) 
330  CONTINUE 

FUN Z =GFUNC * TERMC * C 1 * C 2 + C 3 / GAMMA 

RETURN 

800     FUNZ=GFUNC*TERMC*2 . 0*EXP (GAMMA*Z) /GAMMA 
RETURN 

900     FUNZ=GFUNC* TERMC *EXP  (VS+GAMMA*Z) / (COSH (VS) * GAMMA) 
RETURN 

1000  FUNZ=(  K3/  Kl)*  Ll+(  K3/  K2)*(  L3+  L2+Z) 
RETURN 
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'************************************************ 
:  BOTTOM  LAYER  CALCULATION 

:  THIS   IS  THE  BOTTOM  LAYER  CALCULATION  WHICH  IS  USED  IF  Z  FALLS 

J  INTO  THE  BOTTOM  LAYER 

************************************ 

1500  IF    (GAMMA. EQ . 0.0)   GO  TO  2000 

IF  (VS.GT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0.AND. 
1  VB.GT.5.0)   GO  TO  1900 

IF(VB. GT. 5.0. AND. VS. GT.5.0.AND.VC.LT. 5.0. AND. VI. LT. 5. 0)GO  TO  2100 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2200 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO  TO  2300 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2400 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO  TO  2500 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0)GO  TO  2600 

IF    (VS.LT.5.0)   GO  TO  1550 

Cl=2 . 0*EXP (-VS) 

GO  TO  1560 
1550  C1=1.0/COSH(VS) 
1560  CONTINUE 

IF    (VC.LT.5.0)   GO  TO  157  0 

C2=2 .0*EXP (-VC) 

GO  TO  1580 
1570  C2=1.0/COSH(VC) 
1580  CONTINUE 

IF    (VI.LT.5.0)   GO  TO  1590 

C3=2 . 0*EXP (-VI) 

GO  TO  1600 
1590  C3=1.0/COSH(VI) 
1600  CONTINUE 

C4=SINH (VB) 

FUN Z =GFUNC * C 1 * C 2 * C 3 * C 4 / GAMMA 
RETURN 

1900  FUNZ=GFUNC * 4. 0*EXP (GAMMA* Z) /GAMMA 
RETURN 

2000  FUNZ=(  K3/  Kl)*(  L3+  L2+  Ll+Z) 
RETURN 

2100  FUNZ=GFUNC*EXP (VB-VS) / ( GAMMA* COSH (VC) *COSH(VI) ) 
RETURN 

2200  FUNZ=GFUNC*EXP (VB-VC) / (GAMMA* COSH (VS) *COSH (VI) ) 
RETURN 

2300  FUNZ=GFUNC*EXP (VB-VI) / (GAMMA*COSH (VS) *COSH(VC)  ) 
RETURN 

2400  FUNZ=GFUNC*2 . 0*EXP (VB-VS-VC) / (GAMMA*COSH (VI)  ) 
RETURN 

2500  FUNZ=GFUNC*2 . 0*EXP (VB-VS-VI) / (GAMMA* COSH (VC) ) 
RETURN 

2600  FUNZ=GFUNC*2 . 0*EXP (VB-VC-VI) / (GAMMA* COSH (VS) ) 
RETURN 
END 
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Q*  ****************************************************************************  * 

C  TML  VERSION  1.0  -  VERSION  DATE  06/27/95 
C 

C  THIS  IS  THE  THERMAL  MULTILAYER  PROGRAM  WHICH  IS  BASED  UPON  A  THERMAL 

C  RECURSION  RELATION  WHICH  IS  SIMILAR  TO  THAT  USED  IN  THE  SOLUTION  OF 

C  THE  LAPLACE  EQUATION  FOR  TWO-PROBE  AND  FOUR-PROBE  RESISTANCE  ANALYSIS. 

C  THIS  PROGRAM  CALCULATES  THE  STEADY-STATE  SURFACE  TEMPERATURE  FOR  A 

C  MULTILAYER  RECTANGULAR  STRUCTURE  WITH  AN  ARBITRARY  NUMBER  OF  LAYERS. 

C  THE  CALCULATION  IS  PERFORMED  FOR  AN  ARBITRARY  NUMBER  OF  RECTANGULAR 

C  HEAT  SOURCES/SINKS  ON  THE  TOP  SURFACE. 

C  THE  TEMPERATURE  MAY  BE  CALCULATED  AS  A  POINT  FUNCTION    (X, Y) ,    A  LINE 

C  AVERAGE,    OR  AN  AREA  AVERAGE  ON  THE  TOP  SURFACE. 

C  THE  CALCULATION  FOLLOWS  FROM  THE   INPUT  OF  THE  THICKNESSES  AND  THERMAL 

C  CONDUCTIVITIES  OF  ALL  OF  THE  LAYERS  IN  THE  STRUCTURE. 

C     -       IT  IS  IMPORTANT  TO  EMPHASIZE  THAT  THE  CALCULATION  IS  GENERAL  FOR  THE 

C  MULTILAYER  STRUCTURE  AND  THE  APPLICATION  TO  SEMICONDUCTOR  STRUCTURES 

C  IS  A  SPECIAL  CASE. 
C 

C  THE  STARTING  POINT  IS  GIVEN  IN  EQUATIONS    (13) -(23),    WITH  S=0 

C  (ZERO  FREQUENCY,    STEADY-STATE  CONDITION),    IN  THE  PAPER  BY  KOKKAS . 

C 

C  REFERENCES:    THE  ORIGINAL  MATHEMATICAL  ANALYSIS  OF  THE  THREE-LAYER 

C  STRUCTURE  WAS  PERFORMED  IN  THE  PAPER  "THERMAL  ANALYSIS 

C  OF  MULTIPLE-LAYERED  STRUCTURES"  BY  ACHILLES  G.  KOKKAS, 

C  IEEE  TRANS.   ELEC .   DEV.   VOL.   ED-21,   NO.    11,    674-681  (1974). 

C  THIS  PAPER  WAS  DRAWN  FROM  HIS  PHD  THESIS:   A.   G.  KOKKAS, 

C  "ANALYSIS  AND  DESIGN  OF  ELECTROTHERMAL  INTEGRATED  CIRCUITS, " 

C  PH.D.   THESIS,   MIT,  1972. 

C 

C  THE  ORIGINAL  FORTRAN  IMPLEMENTATION  OF  THE  STEADY  STATE 

C  KOKKAS  EQUATIONS  IS  CONTAINED  IN  THE  TXYZ  CODE  AND 

C  WAS  PRESENTED   IN  THE  REPORT  "SEMICONDUCTOR  MEASUREMENT 

C  TECHNOLOGY:    TXYZ:   A  PROGRAM  FOR  SEMICONDUCTOR  IC  THERMAL 

C  ANALYSIS"     BY  JOHN  ALBERS,    NBS  SPECIAL  PUBLICATION  400-7  6 

C  (APRIL  1984) . 

C 

C  VERSION  2.0  OF  THE  TXYZ  CODE  WAS  PRESENTED   IN  THE  REPORT 

C  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:   VERSION  2 . 0  OF  THE 

C  TXYZ  THERMAL  ANALYSIS  PROGRAM:    TXYZ20"  BY  JOHN  ALBERS, 

C  NIST  SPECIAL  PUBLICATION  400-89    (JUNE  1992) . 
C 

C  VERSION  3.0  OF  THE  TXYZ  CODE  AND  THE  ACCOMPANYING  THERMAL 

C  MULTILAYER  CODES  ARE  DISCUSSED  IN  THIS  REPORT 

C  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:   HOTPAC :  PROGRAMS 

C  FOR  THERMAL  ANALYSIS   INCLUDING  VERSION  3 . 0  OF  THE  TXYZ 

C  PROGRAM,    TXYZ30,   AND  THE  THERMAL  MULTILAYER  PROGRAM,  TML" 

C  BY  JOHN  ALBERS,    NIST  SPECIAL  PUBLICATION  400-96. 

C 

C  THE  REVIEW  AND  APPLICATION  OF  THE  RECURSION  RELATION 

C  TECHNIQUE  FOR  THE  ANALYSIS  OF  THE  LAPLACE  EQUATION  USED  IN 

C  TWO-PROBE  AND  FOUR-PROBE  RESISTANCE  ARE  CONTAINED  IN  THE 

C  REPORT  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:   A  COLLECTION 

C  OF  COMPUTER  PROGRAMS  FOR  TWO-PROBE  RESISTANCE  (SPREADING 

C  RESISTANCE)   AND  FOUR-PROBE  RESISTANCE  CALCULATIONS,  RESPAC 

C  BY  JOHN  ALBERS  AND  HARRY  L.    BERKOWITZ,    NIST  SPECIAL 

C  PUBLICATION  400-91,  1993. 

C 
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C  THE  THERMAL  RECURSION  RELATION  TECHNIQUE  HAS  BEEN  DISCUSSED 

C  IN  THE  PAPERS:   "AN  EXACT  SOLUTION  OF  THE  STEADY-STATE 

C  SURFACE  TEMPERATURE  FOR  A  GENERAL  MULTILAYER  STRUCTURE" 

C  BY  JOHN  ALBERS,    PROCEEDINGS  TENTH  IEEE  SEMI-THERM  SYMPOSIUM, 

C  PP.   129-137,    1994  AND 

C  "AN  EXACT  RECURSION  RELATION  SOLUTION  FOR  THE  STEADY-STATE 

C  SURFACE  TEMPERATURE  OF  A  GENERAL  MULTILAYER  STRUCTURE" 

C  BY  JOHN  ALBERS,    IEEE  TRANSACTIONS  ON  COMPONENTS,  PACKAGING 

C  AND  MANUFACTURING  TECHNOLOGY  -  PART  A,    VOL.    18,   NO.  1, 

C  PP.   31-38,  1995. 

C 

Q*  ****************************************************************************  * 

C  IN  THE  PRESENT  FORM,    THE  PROGRAM  ALLOWS  UP  TO  500  TERMS  TO  BE  INCLUDED 

C  IN  BOTH  THE  N  SUM   (ALONG  X)   AND  THE  M  SUM   (ALONG  Y) .      TO  GO  BEYOND 

C  THIS  NUMBER,    SUBSTITUTE  THE  FOLLOWING  TWO  LINES  WITH  THE  APPROPRIATE 

C  VALUES  OF  NX  AND  MY  FOR  THE  FIRST  TWO  DIMENSION  STATEMENTS  -  MAKE,  SURE  TO 

C  REMOVE  THE  COMMENTS  FROM  THE  NEW  LINES  AND  COMMENT  OUT  THE  REPLACED  LINES. 

C.  DIMENSION  X(100),   Y(100),    Z(100),  COSYT(MY) 

C  DIMENSION  ARUZER (NX, MY) ,    ARFUNZ (NX,  MY) 

C  ALSO  REPLACE  THE  TESTING  LINE  ABOVE  THE  LINE  LABELLED   110  WITH 

C  IF    (NUP.GT.NX.OR.MUP.GT.MY)    GO  TO  3999 


DIMENSION  X (100) ,    Y(100),  COSYT(500) 

DIMENSION  ARUZER(500,  500) ,    ARFUN0  (500,  500) 

DIMENSION  WTSOUR (50) ,    XSOUR(50),  YSOUR(50) 

REAL  LXSOUR(50),    LYSOUR(50),    K(20),    LX,    LY,  L(20) 

DIMENSION  XLINE (30) , XAREA (30) , YAREA (30) 

REAL  LXLINE (30) , LXAREA (30) , LYAREA (30) 

COMMON     LX,    LY,    NLAY,    K,  L 

COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 
PI=3. 14159265 


C  INPUT  DATA  IS  READ  FROM  10  AND  OUTPUT  GOES  TO  12.    THESE  ARE  WRITTEN  IN 

C  LOWER  CASE.     MANY  OPERATING  SYSTEMS  TAKE  UPPER  AND  LOWER  CASE  AS 

C  EQUIVALENT.      THE  UNIX  OPERATING  SYSTEM  VIEWS  THE  UPPER  AND  LOWER 

C  CASE-NAMED  FILES  AS  DIFFERENT  FILES. 

C  UNIX  DEFAULTS  TO  THE  LOWER  CASE  WHICH  IS  USED  HERE. 

open (  unit=10,   file=' input .dat' ,   status  =  'unknown') 
open (  unit=12,   file=' output .dat' ,   status  =  'unknown'  ) 

2  FORMAT (IX,  '  STEADY-STATE  THERMAL  MULTILAYER  CALCULATION'/ 

1'    USING  THE  THERMAL  RECURSION  RELATION  IN  KOKKAS  EQUATIONS') 

3  FORMAT (IX, ' THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  LAYERS') 

4  FORMAT (IX, 'NUMBER  OF  LAYERS  IN  STRUCTURE^ , 12 ) 

5  FORMAT (IX, ' L' , 12, ' =  ',F10.5,'      K' , 12, ' =  ' , F10 . 5) 

6  FORMAT (IX, 'UPPER  SUMMATION  LIMITS  ' , 2X, '   NUP=' , 15, 
1  '  MUP=',I5) 

7  FORMAT (IX, 'NUMBER  OF  HEAT  SOURCES='  , 15) 

8  FORMAT (IX, 'WEIGHTS,    COORDINATES,    LENGTHS,    WIDTHS  OF  SOURCES') 

9  FORMAT ( IX, ' HEAT  SOURCE' , 4X, ' WTSOUR' , 8X, ' XSOUR'  ,  8X,  '  YSOUR'  ,  7X, 
1   ' LXSOUR' , 7X, ' LYSOUR' ) 

10  FORMAT (4X, 13, 5X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5) 

11  FORMAT (IX, 'POWER  DENSITY=' ,F1 1.6) 

27       FORMAT ( IX, ' LX=  ',F7.2,3X, '    LY=  ',F7.2) 

22       FORMAT (IX, F12 . 4, 2X, F12 . 4, 2X, F12 . 4, 2X, F12 .4) 

31  FORMAT (IX, 617) 

32  FORMAT (IX, 'NUMBER  OF  LINES=' , 13) 

33  FORMAT (IX, 13, 3X, F10 . 5,  3X,  F10 . 5) 

34  FORMAT (IX, 'NUMBER  OF  AREAS=' , 13) 
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35  FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X,  F10 . 5) 

40  FORMAT (IX, 'POINT  FUNCTION  EVALUATION  OF  SURFACE  T (X,  Y)  '  ) 

41  FORMAT (IX, ' LINE  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE') 

42  FORMAT ( IX, ' AREA  AVERAGE  EVALUATION  OF  SURFACE  TEMPERATURE') 

43  FORMAT (6X, 'X' , 12X, ' Y' , 11X, ' T (X, Y) ' ) 

44  FORMAT(1X,F10.5, 3X,F10.5, 3X,F10.5) 

45  FORMAT (IX, ' LINE  #' , 5X, ' XLINE' , 7X, ' LXLINE' , 8X, ' Y' ,  9X,  '  AVE  TEMP') 
4  6  FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X,F10.5) 

47  FORMAT (IX, 'AREA  #' , 5X, ' XAREA' , 7X, ' LXAREA' , 8X, ' YAREA' , 7X, 
1   ' L YAREA' , 3X, ' AVE   TEMP' ) 

48  FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, IX,  F10 . 5) 

51  FORMAT (IX, 14, 3X, 14) 

52  FORMAT (IX,  F10 . 5,  3X,F10.5) 

53  FORMAT (IX, II) 

54  FORMAT  (IX, 14, 3X, F10 . 5, 3X,  F10 . 5) 

55  FORMAT (IX, F10 . 5) 

56  FORMAT ( IX, 12 , 3X, F10 . 5 ) 

57  FORMAT(1X,F10.5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5) 

88  FORMAT (IX, ' YOUR  UPPER  LIMIT  OF  SUMMATION  IS  TOO  LARGE.    TRY  AGAIN') 
q* ***************************************************************************** 

C 

C   INPUT  SECTION  

C  THE  FOLLOWING  VARIABLES  ARE  READ  BY  THE  PROGRAM  FROM  FOR010 
C 

C  I  TYPE  TYPE  OF  ANALYSIS 

C  =1  FOR  POINT  FUNCTION,    T(XY, 0),  EVALUATION 

C  =2  FOR  LINE  AVERAGE,    <T(Y,0)>,  EVALUATION 

C  =3  FOR  AREA  AVERAGE,    <T(0)>,  EVALUATION 

C 

C  LX  X  DIMENSION  OF  THE  MULTILAYER  STRUCTURE 

C  LY  Y  DIMENSION  OF  THE  MULTILAYER  STRUCTURE 

C  NLAY  NUMBER  OF  LAYERS   IN  THE  STRUCTURE    (UP  TO  2  0  LAYERS 

C  BUT  THE  CODE  CAN  GO  BEYOND  THIS  BY  EDITING  THE 

C  APPROPRIATE  DIMENSION  STATEMENTS  IN  THE  MAIN  PROGRAM 

C  AND  SUBROUTINES/FUNCTIONS) 

C 

C  THE  NEXT  NLAY  LINES  CONTAIN  THE  THICKNESSES  AND  THERMAL 

C  CONDUCTIVITIES  OF  THE  N-LAYERS  BEGINNING  AT  THE  SURFACE  AND  GOING 

C  TO  THE  BOTTOM  LAYER   (CODE  TAKES  CARE  OF  NUMBERING  FROM  NLAY, ...,1) 

C 

C  L(I),    K (I) -THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  I-TH  LAYER 

C  NUP  UPPER  LIMIT  OF  N  SUM,    X  DIRECTION 

C  MUP  UPPER  LIMIT  OF  M  SUM,    Y  DIRECTION 

C 

C  IEDGEX  INDEX  FOR  GENERATING  THE  VALUES  OF  X  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED  INCREMENT  X  VALUES 

C  =2   TO  READ   IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 

C  IF   IEDGEX=1  THEN  READ  THE  THREE  VARIABLES    (ON  SAME  LINE) 

C  ILX     THE  NUMBER  OF  X  VALUES  TO  BE  USED 

C  XI       THE  VALUE  OF  THE  FIRST  POINT  IN  X 

C  STEPX    (THE   INCREMENT  IN  X) 

C  IF   IEDGEX=2  THEN  READ  THE  VARIABLE  AND  ARRAY    (ONE  PER  LINE) 

C  ILX     THE  NUMBER  OF  X  VALUES  TO  BE  USED 

C  X(I)    THE  ARRAY  OF  X  VALUES    (1=1, ILX) 

C 
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C  IEDGEY  INDEX  FOR  GENERATING  THE  VALUES  OF  Y  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED   INCREMENT  Y  VALUES 

C  =2   TO  READ   IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 

C  IF   IEDGEY=1   THEN  READ  THE  THREE  VARIABLES    (ON  SAME  LINE) 

C  ILY     THE  NUMBER  OF  Y  VALUES  TO  BE  USED 

C  Yl       THE  VALUE  OF  THE  FIRST  POINT  IN  Y 

C  STEPY   (THE  INCREMENT  IN  Y) 

C  IF   IEDGEY=2   THEN  READ  THE  VARIABLE  AND  ARRAY    (ONE  PER  LINE) 

C  ILY     THE  NUMBER  OF  Y  VALUES  TO  BE  USED 

C  Y(I)    THE  ARRAY  OF  Y  VALUES    (1=1, ILY) 
C 

C  NOTE:    THERE   IS  NO  INPUT  FOR  THE  Z  VALUES  AS  Z=0  THROUGHOUT  THE 

C  CALCULATION 

C 

C  IMPORTANT:    THE  POINT  FUNCTION  EVALUATION  OF  THE  TEMPERATURE  IS 

C  THE  MOST  ELEMENTAL  VERSION.      IN  ORDER  TO   SIMPLIFY  THE 

C  CONSTRUCTION  OF  THE  DATA  FILES  FOR  LINE  AND  AREA 

C  AVERAGES,    THE  PROGRAM  EXPECTS  TO  SEE  THE  ABOVE  IEDGEX, 

C  AND   IEDGEY  DATA.      THIS   IS  READ  EVEN   IF   IT   IS  NOT 

C  USED  FOR  THE  AVERAGE  VERSIONS.      HOWEVER,    THE  LINE  AND 

C  AREA  INFORMATION  MAY  THEN  BE  SIMPLY  APPENDED  TO  THE  END 

C  OF  THE  DATA  FILE   IN  ORDER  TO  RUN  THESE  VERSIONS.      SEE  THE 

C  SAMPLE   I/O  FILES  FOR  AN  EXAMPLE  OF  THIS. 

C 

C  NSOUR  NUMBER  OF  HEAT  SOURCES    (UP  TO  50) 

C  P0  POWER  DENSITY 

C 

C  THE  NEXT  NSOUR  LINES  READ  THE  FOLLOWING  INFORMATION  FOR  THE 

C  HEAT  SOURCES    (WITH  ALL  THE   INFORMATION  FOR  EACH  OF  THE  ELEMENTS 

C  ON  A  SINGLE  LINE) 

C 

C  WTSOUR (I) -WEIGHTING  FACTOR  OF   I-TH  SOURCE 

C  (POSITIVE  FOR  SOURCE,    NEGATIVE  FOR  SINK) 

C  XSOUR(I) — X  COORDINATE  OF  ORIGIN  OF  I-TH  SOURCE 

C  LXSOUR (I) -LENGTH  ALONG  X  DIRECTION  OF   I-TH  SOURCE 

C  YSOUR(I) — Y  COORDINATE  OF  ORIGIN  OF  I-TH  SOURCE 

C  LYSOUR (I) -LENGTH  ALONG  Y  DIRECTION  OF   I-TH  SOURCE 

C 

C  IF  ITYPE=1,    THE  POINT  FUNCTION  CALCULATION  CONTINUES  WITH  THE  ABOVE 

C  SET  OF  X,Y  VALUES 

C  IF   ITYPE=2,    THE  LINE  AVERAGE  CALCULATION  READS  THE  FOLLOWING: 

C  NLINE  THE  NUMBER  OF  LINE  SEGMENTS  TO  DO  THE  AVERAGE 

C  THE  NEXT  NLINE  LINES  THEN  READ 

C  XLINE(J),    LXLINE(J) — THE  LOCATION  AND  LENGTH  OF  THE 

C  J-TH  LINE  ELEMENT 

C  IF   ITYPE=3,    THE  AREA  AVERAGE  CALCULATION  READS  THE  FOLLOWING: 

C  NAREA  THE  NUMBER  OF  AREA  SEGMENTS  TO  DO  THE  AVERAGE 

C  THE  NEXT  NAREA  LINES  THEN  READ 

C  XAREA (J) , LXAREA ( J) , YAREA (J) , LYAREA ( J) -THE  LOCATION 

C  AND  LENGTHS  OF  THE  J-TH  AREA  ELMENT 
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C  READ  ITYPE 

READ  (10,*)  ITYPE 

C  READ   LX  AND   LY    (THE  X  AND  Y  DIMENSIONS  OF  THE  RECTANGULAR  STRUCTURE) 

READ  (10,*)    LX,  LY 
C  READ  THE  NUMBER  OF  LAYERS  IN  THE  STRUCTURE 

READ (10,*)  NLAY 

C  READ  LAYER  THICKNESS  AND  THERMAL  CONDUCTIVITY  FROM  SURFACE  TO  BOTTOM 

DO  101   I=NLAY, 1,-1 
READ  (10,  *  )  L(I),K(I) 
101  CONTINUE 

C  READ  NUP  AND  MUP    (UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE   INDEX  N  (X-DIR) 

C  UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE   INDEX  M    (Y-DIR) ) 

C  NUP  AND  MUP  MUST  BE  LESS  THAN  OR  EQUAL  TO  THE  DIMENSIONALITY  OF 

READ  (10,  *) NUP, MUP 

IF    (NUP. GT. 500. OR. MUP. GT. 500)   GO  TO  3999 
READ (10, *) IEDGEX 
GOTO   (110, 115) IEDGEX 

110  READ    (10, *) ILX, XI, STEPX 
DO  111   1=1, ILX 
X(I)=X1+ (1-1) *STEPX 

111  CONTINUE 
GOTO  119 

115  READ  (10,*) ILX 
DO  116  1=1, ILX 
READ (10, *) X  (I) 

116  CONTINUE 

119  READ (10, *) IEDGEY 
GOTO    (120,125)  IEDGEY 

120  READ    (10, *) ILY, Yl, STEPY 
DO  121  1=1, ILY 
Y(I)=Y1+(I-1) *STEPY 

121  CONTINUE 
GOTO  12  9 

125     READ    (10,*) ILY 

DO  126  1=1, ILY 

READ  (10, *) Y(I) 
12  6  CONTINUE 
129  Z=0.0 

C  READ  THE  NUMBER  OF  HEAT  SOURCES  AND  THE  POWER  DENSITY 

C  NOTE-POWER  DENSITY  IS  MULTIPLICATIVE  FACTOR  USUALLY  SET  EQUAL  TO  UNITY 

C  P0   IS  THE  POWER  DENSITY,    ASSUMED  UNIFORM  FOR  ALL  HEATERS 

C  NSOUR  IS  THE  TOTAL  NUMBER  OF  HEATING  ELEMENTS  ON  THE  SURFACE  OF  THE 

C  THE  TOP  LAYER    (UP  TO  50) 

139  READ (10, *)NSOUR, P0 

C  THE  NEXT  LOOP  READS   IN  THE  COORDINATES  OF  THE  ORIGIN  OF  THE 

C  HEATING  ELEMENTS  ALONG  WITH  THEIR  LENGTHS  AND  WIDTHS 

C  THE  WEIGHTING  FACTOR  IS  ALSO  ENTERED    (THIS  IS  REAL,   NON INTEGER) 

C  WTSOUR(I)    IS  THE  WEIGHTING  FACTOR  FOR  THE  I-TH  HEATER  ELEMENT 

C  XSOUR(I)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 

C  LXSOUR(I)    IS  THE  LENGTH  OF  THE   I-TH  HEATER  ALONG  THE  X  DIRECTION 

C  YSOUR(I)    IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 

C  LYSOUR(I)    IS  THE  LENGTH  OF  THE   I-TH  HEATER  ALONG  THE  Y  DIRECTION 

DO   140   1=1, NSOUR 

READ (10,  *)WTSOUR(I) , XSOUR(I) ,    LXSOUR ( I ) , YSOUR ( I ) ,  LYSOUR(I) 

140  CONTINUE 
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C  THIS   IF .. .THEN. . .ELSE   IF  CONSTRUCTION  IS  USED  TO  READ  IN  THE  DATA  FOR 

C  THE  LINES  OR  AREAS  TO  BE  CONSIDERED   IN  THE  CALCULATION 

C  IF  THE  ANALYSIS   IS  FOR  A  POINT  FUNCTION,    THEN  GO  OUT  OF  THE   IF. THEN. ELSE 

IF    (ITYPE.EQ.l)  THEN 
GOTO  17  0 

C  IF  THE  ANALYSIS   IS  FOR  A  LINE  AVERAGE,    THEN  READ  THE  NUMBER  OF  LINE 

C  SEGMENTS  AND  THEIR  LOCATION 

ELSE   IF    ( I TYPE .EQ . 2 )  THEN 
C  READ  THE  NUMBER  OF  LINE  SEGMENTS 

READ (10, * ) NLINE 

DO  150     J=l, NLINE 
C  READ  THE  ORIGIN  AND  LENGTH  OF  EACH  LINE  SEGMENT 

C  NOTE  THAT  THE  AVERAGE   IS  ALONG  THE  X-DIRECTION  FOR  GIVEN  Y  VALUES 

C  TO  DO  THE  AVERAGE  ALONG  THE  Y-DIRECTION  FOR  GIVEN  X,    SIMPLY  ROTATE 

C  THE  STRUCTURE  BY  90  DEGREES  AND  USE  THE  CORRESPONDING  NEW  X' S  AND  Y' S 

C  XLINE(J)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF   J-TH  LINE  ELEMENT 

C  LXLINE(J)    IS  THE  LENGTH  OF  THE  J-TH  LINE  ALONG  THE  X  DIRECTION 

READ (10, * ) XLINE (J) , LXLINE (J) 
150  CONTINUE 

C  IF  THE  ANALYSIS   IS  FOR  AN  AREA  AVERAGE,    THEN  READ  THE  NUMBER  OF  AREAS 

C  AND  THEIR  LOCATIONS 

ELSE   IF    (ITYPE.EQ.3)  THEN 
C  READ  THE  NUMBER  OF  AREAS 

READ (10, * ) NAREA 

DO   160   J=l, NAREA 
C  READ   THE  FOLLOWING 

C  XAREA(J)    IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  J-TH  AREA  ELEMENT 

C  LXAREA(J)    IS  THE  LENGTH  OF  THE  J-TH  AREA  ALONG  THE  X  DIRECTION 

C  YAREA(J)    IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF  J-TH  AREA  ELEMENT 

C  LYAREA(J)    IS  THE  LENGTH  OF  THE  J-TH  AREA  ALONG  THE  Y  DIRECTION 

READ (10, * ) XAREA (J) , LXAREA (J) , YAREA ( J) , LYAREA ( J) 
160  CONTINUE 

END  IF 
17  0  CONTINUE 

£****************************************************************************** 

C  END  OF  DATA   INPUT  SECTION 

C  BEGIN  CALCULATION  OF  T (X, Y, Z) 

C  THE  SUBROUTINES  USED   IN  THE  CALCULATION  ARE: 

C  1)  UZERO(N,M)  -  CALCULATES  THE  FOURIER  COSINE  TRANSFORM  OF  THE 
C  FUNCTION,    U(X,Y),    THE  POWER  DENSITY  FUNCTION  FOR  ALL  OF  THE 

C  HEAT  SOURCES. 

C  2)    FUN0(N,M)    -  CALCULATES  THE  FOURIER  COEFFICIENTS  FOR  THE  LAYERED 
C  STRUCTURE  USING  THE  THERMAL  RECURSION  RELATION 

Q* ***************************************************************************** 

P04LK  =  4.0   *  P0   /    (   LX*  LY) 
PILX  =  PI   /  LX 
PILY  =  PI   /  LY 

C  CALCULATE   THE  FOURIER  COMPONENTS  OF  THE  HEAT   SOURCES,  U(N,M) 

q* ***************************************************************************** 

DO  300  MM=1,MUP 
M  =  MM  -  1 
DO  2  50  NN=1,NUP 
N  =  NN  -  1 

ARUZER (NN, MM) =UZERO (N, M) 
250  CONTINUE 
300  CONTINUE 

£*********************************+++++*+++*+*+++++++++++++*++++*+++++++******* 

C  END  OF  U(N,M)    CALCULATION  AND  BEGINNING  OF  MAJOR  LOOP  FOR  Z 
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q* ********************************************************************* 
C  CALCULATE  THE  FUNO(N,M)    FROM  THE  RECURSION  RELATION 

Q*  ********************************************************************* 

DO  400  MM=1,MUP 
M  =  MM  -  1 
DO  350  NN=1,NUP 
N  =  NN  -  1 

ARFUNO (NN, MM) =FUN0 (N, M) *ARUZER (NN, MM) 
350  CONTINUE 
400  CONTINUE 

C  THE  FOLLOWING  IF. THEN. ELSE  CONSTRUCTION  OPERATES  ACCORDING  TO  THE 

C  TYPE  OF  ANALYSIS  TO  BE  USED. 

C  FOR  THE  POINT  FUNCTION  ANALYSIS,    THE  3000  LOOP   IS  USED 

C  FOR  THE  LINE  AVERAGE  ANALYSIS,    THE  4000  LOOP  IS  USED 

C  FOR  THE  AREA  AVERAGE  ANALYSIS,    THE  5000  LOOP   IS  USED 

C  THIS  PORTION  DOES  THE  POINT  FUNCTION  CALCULATION 

IF    (ITYPE.EQ.l)  THEN 
C  BEGINNING  OF  POINT  FUNCTION  ANALYSIS 

WRITE  (12, 2) 

WRITE  (12, 40) 

WRITE  (12,  43) 

DO  3000     IY=1, ILY 

DO  3100  MM=1,MUP 
M  =  MM  -  1 

COSYT (MM) =COS (FLOAT (M) *Y(IY) *PILY) 
3100  CONTINUE 

DO  3000   IX=1, ILX 
SUM=0 . 0 

DO  3300  MM=1,MUP 

M  =  MM  -  1 

DO  3200  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF  (N.EQ.0)  NDN=1 
IF    (M.EQ.0)  NDM=1 

TOP  =  ARFUNO (NN, MM)    *  COS  (FLOAT (N) *X ( IX) *PILX)    *  COSYT (MM) 
BOTTOM= (NDN+1) * (NDM+1) 
TSUM=TOP /BOTTOM 
SUM=SUM+TSUM 
32  00  CONTINUE 
3300  CONTINUE 

TEMP  =  P04LK  *  SUM/K (NLAY) 
WRITE (12, 44)X(IX),Y(IY), TEMP 
3000  CONTINUE 
C  END  OF  POINT  FUNCTION  CALCULATION 

ELSE  IF    (I TYPE . EQ . 2 )  THEN 
C  THIS  PORTION  DOES  THE  LINE  AVERAGE  CALCULATION 

WRITE  (12, 2) 
WRITE (12, 41) 
WRITE (12, 45) 
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DO  4000  IY=1,ILY 
DO  4010  MM=1,MUP 
M  =  MM  -  1 

COSYT (MM) =COS (FLOAT (M) *Y(IY) *PILY) 
4010  CONTINUE 

DO  4000  J=1,NLINE 
SUM=0 . 0 

DO  4200  MM=1,MUP 

M  =  MM  -  1 

DO  4100  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF    (N.EQ.0)  NDN=1 
IF    (M.EQ.0)  NDM=1 
IF(N.EQ.O)   GO  TO  4160 
TERMX= SIN (FLOAT (N) *PI* (XLINE (J) +LXLINE (J)  )  /LX) 
1  -SIN(FLOAT(N)*PI*XLINE(J)/LX) 
TERMX=TERMX*LX/ (FLOAT (N) *PI) 
GO  TO  4165 
4160     TERMX=LXLINE (J) 
4165  CONTINUE 

TOP  =  ARFUN0 (NN, MM)    *   COSYT (MM)    *  TERMX 
BOTTOM= (NDN+1) * (NDM+1) 
TSUM=TOP /BOTTOM 
SUM=SUM+TSUM 
4100  CONTINUE 

TEMP  =  P04LK  *   SUM  /    (LXLINE (J)    *  K(NLAY)) 
42  00  CONTINUE 

WRITE (12, 46) J, XLINE (J) , LXLINE (J) , Y(IY) , TEMP 
4  000  CONTINUE 
C  END  OF  LINE  AVERAGE  PORTION  OF  THE  CODE 

ELSE  IF    (ITYPE.EQ.3)  THEN 
C  THIS  PART  PERFORMS  THE  AREA  AVERAGE  CALCULATION 

WRITE (12, 2) 
WRITE  (12, 42) 
WRITE  (12, 47) 
DO  5000   J= 1 , NARE A 
SUM=0.0 

DO  5100  MM=1,MUP 

M  =  MM  -  1 

DO  5200  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF    (N.EQ.0)  NDN=1 
IF    (M.EQ.0)  NDM=1 
AREA=0 . 0 

IF (N.EQ.0)   GO  TO  5160 

TERMX  =  SIN  (FLOAT  (N) *PI* (XAREA (J)    +     LXAREA ( J) ) /  LX) 
1  -  SIN (FLOAT (N) *PI*XAREA (J) /  LX) 

TERMX=TERMX *  LX/ (FLOAT (N) *PI) 
GO  TO  5165 
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5160  TERMX=  LXAREA(J) 

5165  IF(M.EQ.O)   GO  TO  5164 

TERMY  =  SIN (FLOAT  (M) *PI*  (YAREA (J)    +     LYAREA ( J) ) /  LY) 
1  -  SIN (FLOAT (M) *P I* YAREA (J) /  LY) 

TERMY=TERMY *  LY/ (FLOAT (M) *P I ) 
GO  TO  5166 
5164  TERMY=  LYAREA (J) 

5166  TERMI=TERMX*TERMY 
AREA=TERMI 

TOP  =  ARFUNO  (NN,  MM)    *  AREA 
BOTTOM= (NDN+1) * (NDM+1) 
TSUM=TOP /BOTTOM 
SUM=SUM+TSUM 
5200  CONTINUE 

TEMP  =  P04LK  *   SUM  /    (LXAREA (J) *LYAREA (J) *K (NLAY) ) 
5100  CONTINUE 

WRITE (12, 48) J, XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA (J) , TEMP 
5000  CONTINUE 
END  IF 

WRITE (12,27)    LX,  LY 
WRITE (12, 4)  NLAY 
WRITE (12, 3) 
DO  3001   I=NLAY, 1,-1 
WRITE  (12,  5)    I,L(I),  I,K(I) 
3001  CONTINUE 

WRITE (12, 6)NUP,MUP 
WRITE (12, 11)  P0 
WRITE (12, 7)NSOUR 
WRITE (12, 8) 
WRITE  (12,  9) 
DO  3888  I=l,NSOUR 

WRITE (12,  10) I, WTSOUR(I) , XSOUR(l) , YSOUR(I) ,    LXSOUR(I),  LYSOUR(I) 
3888  CONTINUE 

IF (ITYPE.EQ. 1)  THEN 

GOTO  7000 
ELSE   IF    (ITYPE.EQ. 2)  THEN 

WRITE (12, 32)NLINE 

DO  6100  J=1,NLINE 

WRITE (12, 33) J, XLINE (J) , LXLINE (J) 
6100  CONTINUE 

ELSE  IF    (ITYPE.EQ. 3)  THEN 
WRITE (12, 34)  NAREA 
DO  6200  J=l, NAREA 

WRITE (12,  35) J, XAREA (J) , LXAREA (J) , YAREA (J)  ,  LYAREA (J) 
6200  CONTINUE 

END  IF 

GO  TO  7000 
3999  WRITE (12,  88) 
7000  STOP 

END 

^*********************************^ 
;  END  OF  THE  MAIN  PROGRAM 
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:  FUNCTION  UZERO  LISTING 

DESCRIPTION  OF  THE  FUNCTION  UZERO (N,M) 

THIS  FUNCTION  CALCULATES  THE  DOUBLE  FOURIER  COSINE  TRANSFORM 
OF  THE  POWER  DENSITY  FUNCTION,    U(X,Y).      THIS  IS  THE  TRANSFORM 
FOR  ALL  OF  THE  HEAT  SOURCES.      THE  ASSUMPTION  IS  MADE  THAT  THE 
POWER  DENSITY  IS  UNIFORM  AND  EQUAL  TO  UNITY  OVER  THE  SURFACE 
OF  THE  HEATING  ELEMENTS.      THAT  IS, 

U(X,Y)=1      (XSOUR(I) <=X<=XSOUR(I) +LXSOUR(I)  AND 

YSOUR ( I ) <=Y<=YXOUR ( I ) +LYSOUR ( I )    ) . 
U(X,Y)=0  OTHERWISE. 
UNDER  THESE  CONDITIONS,    IT  IS  POSSIBLE  TO  ANALYTICALLY  EVALUATE 
THE  DOUBLE   INTEGRAL  FOR  EACH  HEATING  ELEMENT.     AS  THE  HEATING 
ELEMENTS  ARE  ASSUMED  TO  BE   INDEPENDENT,    THE  CONTRIBUTION  FROM 
EACH  ELEMENT  MAY  BE  ADDED  TO  OBTAIN  THE  U (N, M)    FOR  ALL. 


C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


100 
150 


200 
250 

500 


NONUNIFORM  POWER  DENSITIES  MAY  BE  TAKEN  CARE  OF  BY  USING  THE 

WEIGHTING  FACTOR  FOR  EACH  ELEMENT  (READ  IN  THE  MAIN  PROGRAM) 
***************+************★******** 

FUNCTION  UZERO (N,M) 

DIMENSION  WTSOUR( 50) ,    XSOUR(50),  YSOUR(50) 

REAL  LXSOUR(50),    LYSOUR(50),    K(20),    LX,    LY,  L(20) 

COMMON     LX,    LY,    NLAY,    K,  L 

COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 

PI=3 . 14159265 

UZERO=0 . 0 

DO  500  1=1, NSOUR 

IF(N.EQ.O)   GO  TO  100 

TERMX  =  SIN (FLOAT (N) *PI* (XSOUR (I)    +     LXSOUR(I))/  LX) 
1  -  SIN (FLOAT (N) *PI*XSOUR ( I ) /  LX) 

TERMX=TERMX *   LX/  (FLOAT (N) *PI ) 
GO  TO  150 
TERMX=  LXSOUR (I) 
IF(M.EQ.O)   GO  TO  200 

TERMY  =  SIN (FLOAT (M) *PI* (YSOUR (I)    +     LYSOUR(I))/  LY) 
1  -  SIN (FLOAT (M) *P I* YSOUR (I)  /  LY) 

TERMY=TERMY *   LY/ (FLOAT (M) *PI) 
GO  TO  250 
TERMY=  LYSOUR (I) 
TERMI=TERMX* TERMY 
UZERO=UZERO+TERMI *WTSOUR ( I ) 
CONTINUE 
RETURN 
END 

********************************************* 
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C  FUNCTION  FUNO  LISTING 

C 

C  DESCRIPTION  OF  THE  FUNCTION  FUN0(N,M) 

C 

Q* ***************************************************************************** 

FUNCTION  FUN0(N,M) 

DIMENSION  WTSOUR(50),   XSOUR(50),  YSOUR(50) 

REAL  LXSOUR(50),    LYSOUR(50),    K(20),    LX,    LY,  L(20) 

COMMON     LX,    LY,    NLAY,    K,  L 

COMMON  NSOUR, WTSOUR, XSOUR, YSOUR,    LXSOUR,  LYSOUR 
PI=3. 14159265 

GAMMA=SQRT ( (FLOAT  (N) *PI/  LX) **2  +    (FLOAT (M) *PI/  LY)**2) 

IF    (GAMMA. EQ.O)   GO  TO  120 
C  FOR  GAMMA. GT.O,    START  WITH  BOTTOM  LAYER 

FUN0=TANH (GAMMA*L (1) ) 

IF    (NLAY.EQ.l)    GO  TO  110 
C  THE   100  LOOP  DOES  THE  RECURSION  RELATION  FOR  NONZERO  GAMMA 

DO  100  1=2, NLAY 

TOP=FUN0*K(I) +K(I-1) *TANH (GAMMA*L (I) ) 

BOT=K (I— 1) +K (I) *FUN0*TANH (GAMMA*L (I) ) 

FUN0=TOP/BOT 
100  CONTINUE 
110     FUN 0=FUN0 /GAMMA 

RETURN 
120  FUN0=L(1) 

IF    (NLAY.EQ.l)   GO  TO  130 
C  THE  125  LOOP  DOES  THE  RECURSION  RELATION  FOR  GAMMA  EQUAL  TO  ZERO 

DO  125  1=2, NLAY 

FUN0= (FUN0*K(I)+K(I-1) *L(I) ) /K(I-l) 
125  CONTINUE 
130  RETURN 

END 
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Periodical 


Journal  of  Research  of  the  National  Institute  of  Standards  and  Technology — Reports  NIST  research 
and  development  in  those  disciplines  of  the  physical  and  engineering  sciences  in  which  the  Institute  is 
active.  These  include  physics,  chemistry,  engineering,  mathematics,  and  computer  sciences.  Papers  cover  a 
broad  range  of  subjects,  with  major  emphasis  on  measurement  methodology  and  the  basic  technology 
underlying  standardization.  Also  included  from  time  to  time  are  survey  articles  on  topics  closely  related  to 
the  Institute's  technical  and  scientific  programs.  Issued  six  times  a  year. 

Nonperiodicals 


Monographs — Major  contributions  to  the  technical  literature  on  various  subjects  related  to  the 
Institute's  scientific  and  technical  activities. 

Handbooks — Recommended  codes  of  engineering  and  industrial  practice  (including  safety  codes)  devel  - 
oped in  cooperation  with  interested  industries,  professional  organizations,  and  regulatory  bodies. 
Special  Publications — Include  proceedings  of  conferences  sponsored  by  NIST,  NIST  annual  reports,  and 
other  special  publications  appropriate  to  this  grouping  such  as  wall  charts,  pocket  cards,  and  bibliographies. 

National  Standard  Reference  Data  Series — Provides  quantitative  data  on  the  physical  and  chemical 
properties  of  materials,  compiled  from  the  world's  literature  and  critically  evaluated.  Developed  under  a 
worldwide  program  coordinated  by  NIST  under  the  authority  of  the  National  Standard  Data  Act  (Public 
Law  90-396).  NOTE:  The  Journal  of  Physical  and  Chemical  Reference  Data  (JPCRD)  is  published 
bimonthly  for  NIST  by  the  American  Chemical  Society  (ACS)  and  the  American  Institute  of  Physics  (AIP). 
Subscriptions,  reprints,  and  supplements  are  available  from  ACS,  1155  Sixteenth  St.,  NW,  Washington,  DC 
20056. 

Building  Science  Series — Disseminates  technical  information  developed  at  the  Institute  on  building 
materials,  components,  systems,  and  whole  structures.  The  series  presents  research  results,  test  methods,  and 
performance  criteria  related  to  the  structural  and  environmental  functions  and  the  durability  and  safety 
characteristics  of  building  elements  and  systems. 

Technical  Notes — Studies  or  reports  which  are  complete  in  themselves  but  restrictive  in  their  treatment  of 
a  subject.  Analogous  to  monographs  but  not  so  comprehensive  in  scope  or  definitive  in  treatment  of  the 
subject  area.  Often  serve  as  a  vehicle  for  final  reports  of  work  performed  at  NIST  under  the  sponsorship  of 
other  government  agencies. 

Voluntary  Product  Standards — Developed  under  procedures  published  by  the  Department  of  Commerce 
in  Part  10,  Title  15,  of  the  Code  of  Federal  Regulations.  The  standards  establish  nationally  recognized 
requirements  for  products,  and  provide  all  concerned  interests  with  a  basis  for  common  understanding  of 
the  characteristics  of  the  products.  NIST  administers  this  program  in  support  of  the  efforts  of  private-sector 
standardizing  organizations. 

Order  the  following  NIST  publications — FIPS  and  NISTIRs—from  the  National  Technical  Information 
Service,  Springfield,  VA  22 1 61. 

Federal  Information  Processing  Standards  Publications  (FIPS  PUB) — Publications  in  this  series 
collectively  constitute  the  Federal  Information  Processing  Standards  Register.  The  Register  serves  as  the 
official  source  of  information  in  the  Federal  Government  regarding  standards  issued  by  NIST  pursuant  to 
the  Federal  Property  and  Administrative  Services  Act  of  1949  as  amended,  Public  Law  89-306  (79  Stat. 
1127),  and  as  implemented  by  Executive  Order  11717  (38  FR  12315,  dated  May  11,  1973)  and  Part  6  of 
Title  15  CFR  (Code  of  Federal  Regulations). 

NIST  Interagency  Reports  (NISTIR) — A  special  series  of  interim  or  final  reports  on  work  performed  by 
NIST  for  outside  sponsors  (both  government  and  nongovernment).  In  general,  initial  distribution  is  handled 
by  the  sponsor;  public  distribution  is  by  the  National  Technical  Information  Service,  Springfield,  VA  22161, 
in  paper  copy  or  microfiche  form. 
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